programing tip

expr에서 오버플로를 피하는 방법

itbloger 2020. 6. 2. 08:22
반응형

expr에서 오버플로를 피하는 방법 A * B-C * D


다음과 같은 표현식을 계산해야합니다. A*B - C*D, 유형은 다음과 같습니다 . signed long long int A, B, C, D;각 숫자는 실제로 유형이 넘칠 수 없습니다. A*B오버플로가 발생할 수 있지만 동시에 표현 A*B - C*D은 매우 작을 수 있습니다. 올바르게 계산하려면 어떻게해야합니까?

예 : MAX * MAX - (MAX - 1) * (MAX + 1) == 1, where MAX = LLONG_MAX - n및 n-자연수


이것은 너무 사소한 것 같아요. 그러나 A*B넘칠 수있는 것입니다.

정밀도를 잃지 않고 다음을 수행 할 수 있습니다

A*B - C*D = A(D+E) - (A+F)D
          = AD + AE - AD - DF
          = AE - DF
             ^smaller quantities E & F

E = B - D (hence, far smaller than B)
F = C - A (hence, far smaller than C)

이 분해는 추가수행 할 수 있습니다 .
@Gian이 지적했듯이 유형이 오랫동안 부호가 없으면 빼기 연산 중에주의를 기울여야 할 수도 있습니다.


예를 들어, 질문에 해당하는 경우 한 번만 반복하면됩니다.

 MAX * MAX - (MAX - 1) * (MAX + 1)
  A     B       C           D

E = B - D = -1
F = C - A = -1

AE - DF = {MAX * -1} - {(MAX + 1) * -1} = -MAX + MAX + 1 = 1

가장 단순하고 가장 일반적인 해결책은 긴 정수 라이브러리 (예 : http://gmplib.org/ )를 사용하거나 구조체 또는 배열을 사용하여 표현하고 일종의 긴 곱셈을 구현 하여 오버플로 할 수없는 표현을 사용하는 것입니다 ( 즉, 각 숫자를 두 개의 32 비트 반으로 분리하고 다음과 같이 곱셈을 수행합니다.

(R1 + R2 * 2^32 + R3 * 2^64 + R4 * 2^96) = R = A*B = (A1 + A2 * 2^32) * (B1 + B2 * 2^32) 
R1 = (A1*B1) % 2^32
R2 = ((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) % 2^32
R3 = (((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) / 2^32 + (A1*B2) / 2^32 + (A2*B1) / 2^32 + (A2*B2) % 2^32) %2^32
R4 = ((((A1*B1) / 2^32 + (A1*B2) % 2^32 + (A2*B1) % 2^32) / 2^32 + (A1*B2) / 2^32 + (A2*B1) / 2^32 + (A2*B2) % 2^32) / 2^32) + (A2*B2) / 2^32

최종 결과가 64 비트에 적합하다고 가정하면 실제로 대부분의 R3 비트는 필요하지 않으며 R4는 필요하지 않습니다


이는 랩 어라운드 부호있는 오버 플로우에 의존하기 때문에 표준이 아닙니다. (GCC에는이를 가능하게하는 컴파일러 플래그가 있습니다.)

그러나의 모든 계산 long long을 수행하는 경우 수식을 직접 적용한
(A * B - C * D)결과는 올바른 결과가에 맞는 한 정확합니다 long long.


부호없는 정수를 부호있는 정수로 캐스팅하는 구현 정의 동작에만 의존하는 해결 방법이 있습니다. 그러나 이것은 오늘날 거의 모든 시스템에서 작동 할 것으로 예상됩니다.

(long long)((unsigned long long)A * B - (unsigned long long)C * D)

이는 unsigned long long표준에 의해 오버 플로우 동작이 보장되는 곳으로 입력을 캐스트합니다 . 마지막에 부호있는 정수로 다시 캐스팅하는 것은 구현 정의 부분이지만 오늘날 거의 모든 환경에서 작동합니다.


더 많은 pedantic 솔루션이 필요하다면 "긴 산술"을 사용해야한다고 생각합니다


이것은 작동해야합니다 (제 생각에).

signed long long int a = 0x7ffffffffffffffd;
signed long long int b = 0x7ffffffffffffffd;
signed long long int c = 0x7ffffffffffffffc;
signed long long int d = 0x7ffffffffffffffe;
signed long long int bd = b / d;
signed long long int bdmod = b % d;
signed long long int ca = c / a;
signed long long int camod = c % a;
signed long long int x = (bd - ca) * a * d - (camod * d - bdmod * a);

내 파생물은 다음과 같습니다.

x = a * b - c * d
x / (a * d) = (a * b - c * d) / (a * d)
x / (a * d) = b / d - c / a

now, the integer/mod stuff:
x / (a * d) = (b / d + ( b % d ) / d) - (c / a + ( c % a ) / a )
x / (a * d) = (b / d - c / a) - ( ( c % a ) / a - ( b % d ) / d)
x = (b / d - c / a) * a * d - ( ( c % a ) * d - ( b % d ) * a)

모든 값에 대해 가장 큰 공통 요소를 계산 한 다음 산술 연산을 수행 한 다음 다시 곱하기 전에 해당 요소로 나누어 계산할 수 있습니다. 이 같은 요인이 있지만, 있다고 가정합니다 (경우, 예를 들어 A, B, CD프라임은 비교적 될 일이, 그들은 공통의 요소가되지 않습니다).

마찬가지로 로그 스케일 작업을 고려할 수 있지만 수치 정밀도에 따라 약간 무섭습니다.


결과가 long int int 인 경우 A * BC * D 표현식은 산술 모드 2 ^ 64를 수행하므로 괜찮으며 올바른 결과를 제공합니다. 문제는 결과가 긴 long int에 맞는지 아는 것입니다. 이것을 감지하려면 double을 사용하여 다음 트릭을 사용할 수 있습니다.

if( abs( (double)A*B - (double)C*D ) > MAX_LLONG ) 
    Overflow
else 
    return A*B-C*D;

이 접근 방식의 문제점은 복식의 가수 정밀도 (54 비트?)에 의해 제한되므로 제품 A * B 및 C * D를 63 + 54 비트 (또는 약간 더 적음)로 제한해야한다는 것입니다.


E = max(A,B,C,D)
A1 = A -E;
B1 = B -E;
C1 = C -E;
D1 = D -E;

그때

A*B - C*D = (A1+E)*(B1+E)-(C1+E)(D1+E) = (A1+B1-C1-D1)*E + A1*B1 -C1*D1

배열에 각 숫자를 쓸 수 있으며 각 요소는 숫자이며 다항식으로 계산을 수행 할 수 있습니다 . 결과 다항식 (배열)을 가져 와서 배열의 각 요소에 10을 곱하여 배열의 위치 거듭 제곱 (첫 번째 위치가 가장 크고 마지막은 0)으로 결과를 계산합니다.

숫자 123는 다음과 같이 표현 될 수 있습니다.

123 = 100 * 1 + 10 * 2 + 3

이것에 대한 배열을 만듭니다 [1 2 3].

모든 숫자 A, B, C 및 D에 대해이 작업을 수행 한 다음 다항식으로 곱합니다. 결과 다항식이 있으면 그 숫자를 재구성하십시오.


signed long long int의지가 유지되지 않는 동안 A*B그들 중 두 명이됩니다. 따라서 A*B다른 지수의 트리 항으로 분해 될 수 있습니다 signed long long int.

A1=A>>32;
A0=A & 0xffffffff;
B1=B>>32;
B0=B & 0xffffffff;

AB_0=A0*B0;
AB_1=A0*B1+A1*B0;
AB_2=A1*B1;

동일합니다 C*D.

직선 방식 Folowing의 subraction는 모든 쌍 수행 할 수 AB_iCD_i각각에 대한 추가 캐리 비트 (정확하게 1 비트 정수)를 사용하여, 마찬가지로. 따라서 E = A * BC * D라고하면 다음과 같은 결과가 나타납니다.

E_00=AB_0-CD_0 
E_01=(AB_0 > CD_0) == (AB_0 - CD_0 < 0) ? 0 : 1  // carry bit if overflow
E_10=AB_1-CD_1 
...

우리는 상반신 E_10E_20(32만큼 쉬프트하고 추가 한 다음, 상반신을 지우십시오)를 계속해서 전송합니다 E_10.

이제 캐리 비트 E_11를 비 캐리 부분에서 가져온 오른쪽 부호와 함께 추가 하여 캐리 비트 제거 할 수 있습니다 E_20. 이것이 오버플로를 유발하면 결과도 맞지 않습니다.

E_10이제 상위 절반을 E_00(shift, add, erase) 및 carry bit 에서 가져 오기에 충분한 '공간'이 있습니다 E_01.

E_10이제 더 커질 수 있으므로로 전송을 반복합니다 E_20.

이 시점에서 E_200 되어야합니다. 그렇지 않으면 결과가 맞지 않습니다. E_10전송 결과로 인해 상단 도 비어 있습니다.

마지막 단계는 아래쪽 절반 E_20E_10다시 옮깁니다.

보류에 E=A*B+C*D적합한 기대가 있다면signed long long int

E_20=0
E_10=0
E_00=E

최종 결과가 정수 유형으로 표현 가능하다는 것을 알고 있다면 아래 코드를 사용하여이 계산을 빠르게 수행 할 수 있습니다. C 표준은 부호없는 산술이 모듈로 산술이고 오버플로되지 않도록 지정하기 때문에 부호없는 유형을 사용하여 계산을 수행 할 수 있습니다.

다음 코드는 동일한 너비의 부호없는 유형이 있고 부호있는 유형은 모든 비트 패턴을 사용하여 값을 표시한다고 가정합니다 (트랩 표현 없음, 부호있는 유형의 최소값은 부호없는 유형의 계수의 절반에 해당함). 이것이 C 구현에서 유지되지 않는 경우,이를 위해 ConvertToSigned 루틴을 간단하게 조정할 수 있습니다.

다음은 코드를 사용 signed char하고 unsigned char보여줍니다. 귀하의 구현을 위해,의 정의 변경 Signedtypedef signed long long int Signed;와의 정의 Unsigned에를 typedef unsigned long long int Unsigned;.

#include <limits.h>
#include <stdio.h>
#include <stdlib.h>


//  Define the signed and unsigned types we wish to use.
typedef signed char   Signed;
typedef unsigned char Unsigned;

//  uHalfModulus is half the modulus of the unsigned type.
static const Unsigned uHalfModulus = UCHAR_MAX/2+1;

//  sHalfModulus is the negation of half the modulus of the unsigned type.
static const Signed   sHalfModulus = -1 - (Signed) (UCHAR_MAX/2);


/*  Map the unsigned value to the signed value that is the same modulo the
    modulus of the unsigned type.  If the input x maps to a positive value, we
    simply return x.  If it maps to a negative value, we return x minus the
    modulus of the unsigned type.

    In most C implementations, this routine could simply be "return x;".
    However, this version uses several steps to convert x to a negative value
    so that overflow is avoided.
*/
static Signed ConvertToSigned(Unsigned x)
{
    /*  If x is representable in the signed type, return it.  (In some
        implementations, 
    */
    if (x < uHalfModulus)
        return x;

    /*  Otherwise, return x minus the modulus of the unsigned type, taking
        care not to overflow the signed type.
    */
    return (Signed) (x - uHalfModulus) - sHalfModulus;
}


/*  Calculate A*B - C*D given that the result is representable as a Signed
    value.
*/
static signed char Calculate(Signed A, Signed B, Signed C, Signed D)
{
    /*  Map signed values to unsigned values.  Positive values are unaltered.
        Negative values have the modulus of the unsigned type added.  Because
        we do modulo arithmetic below, adding the modulus does not change the
        final result.
    */
    Unsigned a = A;
    Unsigned b = B;
    Unsigned c = C;
    Unsigned d = D;

    //  Calculate with modulo arithmetic.
    Unsigned t = a*b - c*d;

    //  Map the unsigned value to the corresponding signed value.
    return ConvertToSigned(t);
}


int main()
{
    //  Test every combination of inputs for signed char.
    for (int A = SCHAR_MIN; A <= SCHAR_MAX; ++A)
    for (int B = SCHAR_MIN; B <= SCHAR_MAX; ++B)
    for (int C = SCHAR_MIN; C <= SCHAR_MAX; ++C)
    for (int D = SCHAR_MIN; D <= SCHAR_MAX; ++D)
    {
        //  Use int to calculate the expected result.
        int t0 = A*B - C*D;

        //  If the result is not representable in signed char, skip this case.
        if (t0 < SCHAR_MIN || SCHAR_MAX < t0)
            continue;

        //  Calculate the result with the sample code.
        int t1 = Calculate(A, B, C, D);

        //  Test the result for errors.
        if (t0 != t1)
        {
            printf("%d*%d - %d*%d = %d, but %d was returned.\n",
                A, B, C, D, t0, t1);
            exit(EXIT_FAILURE);
        }
    }
    return 0;
}

오버플로하지 않는 더 작은 성분으로 방정식을 깰 수 있습니다.

AB - CD
= [ A(B - N) - C( D - M )] + [AN - CM]

= ( AK - CJ ) + ( AN - CM)

    where K = B - N
          J = D - M

구성 요소가 여전히 오버 플로우되면 작은 구성 요소를 재귀 적으로 분할 한 다음 다시 결합 할 수 있습니다.


모든 경우를 다루지 않았거나 엄격하게 테스트하지는 않았지만 16 비트 CPU에서 32 비트 정수 수학을 시도 할 때 80 년대에 사용했던 기술을 구현합니다. 기본적으로 32 비트를 2 개의 16 비트 장치로 분할하고 별도로 작업합니다.

public class DoubleMaths {
  private static class SplitLong {
    // High half (or integral part).
    private final long h;
    // Low half.
    private final long l;
    // Split.
    private static final int SPLIT = (Long.SIZE / 2);

    // Make from an existing pair.
    private SplitLong(long h, long l) {
      // Let l overflow into h.
      this.h = h + (l >> SPLIT);
      this.l = l % (1l << SPLIT);
    }

    public SplitLong(long v) {
      h = v >> SPLIT;
      l = v % (1l << SPLIT);
    }

    public long longValue() {
      return (h << SPLIT) + l;
    }

    public SplitLong add ( SplitLong b ) {
      // TODO: Check for overflow.
      return new SplitLong ( longValue() + b.longValue() );
    }

    public SplitLong sub ( SplitLong b ) {
      // TODO: Check for overflow.
      return new SplitLong ( longValue() - b.longValue() );
    }

    public SplitLong mul ( SplitLong b ) {
      /*
       * e.g. 10 * 15 = 150
       * 
       * Divide 10 and 15 by 5
       * 
       * 2 * 3 = 5
       * 
       * Must therefore multiply up by 5 * 5 = 25
       * 
       * 5 * 25 = 150
       */
      long lbl = l * b.l;
      long hbh = h * b.h;
      long lbh = l * b.h;
      long hbl = h * b.l;
      return new SplitLong ( lbh + hbl, lbl + hbh );
    }

    @Override
    public String toString () {
      return Long.toHexString(h)+"|"+Long.toHexString(l);
    }
  }

  // I'll use long and int but this can apply just as easily to long-long and long.
  // The aim is to calculate A*B - C*D without overflow.
  static final long A = Long.MAX_VALUE;
  static final long B = Long.MAX_VALUE - 1;
  static final long C = Long.MAX_VALUE;
  static final long D = Long.MAX_VALUE - 2;

  public static void main(String[] args) throws InterruptedException {
    // First do it with BigIntegers to get what the result should be.
    BigInteger a = BigInteger.valueOf(A);
    BigInteger b = BigInteger.valueOf(B);
    BigInteger c = BigInteger.valueOf(C);
    BigInteger d = BigInteger.valueOf(D);
    BigInteger answer = a.multiply(b).subtract(c.multiply(d));
    System.out.println("A*B - C*D = "+answer+" = "+answer.toString(16));

    // Make one and test its integrity.
    SplitLong sla = new SplitLong(A);
    System.out.println("A="+Long.toHexString(A)+" ("+sla.toString()+") = "+Long.toHexString(sla.longValue()));

    // Start small.
    SplitLong sl10 = new SplitLong(10);
    SplitLong sl15 = new SplitLong(15);
    SplitLong sl150 = sl10.mul(sl15);
    System.out.println("10="+sl10.longValue()+"("+sl10.toString()+") * 15="+sl15.longValue()+"("+sl15.toString()+") = "+sl150.longValue() + " ("+sl150.toString()+")");

    // The real thing.
    SplitLong slb = new SplitLong(B);
    SplitLong slc = new SplitLong(C);
    SplitLong sld = new SplitLong(D);
    System.out.println("B="+Long.toHexString(B)+" ("+slb.toString()+") = "+Long.toHexString(slb.longValue()));
    System.out.println("C="+Long.toHexString(C)+" ("+slc.toString()+") = "+Long.toHexString(slc.longValue()));
    System.out.println("D="+Long.toHexString(D)+" ("+sld.toString()+") = "+Long.toHexString(sld.longValue()));
    SplitLong sanswer = sla.mul(slb).sub(slc.mul(sld));
    System.out.println("A*B - C*D = "+sanswer+" = "+sanswer.longValue());

  }

}

인쇄물:

A*B - C*D = 9223372036854775807 = 7fffffffffffffff
A=7fffffffffffffff (7fffffff|ffffffff) = 7fffffffffffffff
10=10(0|a) * 15=15(0|f) = 150 (0|96)
B=7ffffffffffffffe (7fffffff|fffffffe) = 7ffffffffffffffe
C=7fffffffffffffff (7fffffff|ffffffff) = 7fffffffffffffff
D=7ffffffffffffffd (7fffffff|fffffffd) = 7ffffffffffffffd
A*B - C*D = 7fffffff|ffffffff = 9223372036854775807

작동하는 것처럼 보입니다.

나는 부호 넘침 등을 보는 것과 같은 미묘한 부분을 놓친 것이지만 내 생각에는 본질이 있다고 생각합니다.


For the sake of completeness, since no one mentioned it, some compilers (e.g. GCC) actually provide you with a 128 bit integer nowadays.

Thus an easy solution could be:

(long long)((__int128)A * B - (__int128)C * D)

AB-CD = (AB-CD) * AC / AC = (B/C-D/A)*A*C. Neither B/C nor D/A can overflow, so calculate (B/C-D/A) first. Since the final result won't overflow according to your definition, you can safely perform the remaining multiplications and calculate (B/C-D/A)*A*C which is the required result.

Note, if your input can be extremely small as well, the B/C or D/A can overflow. If it's possible, more complex manipulations might be required according to input inspection.


Choose K = a big number (eg. K = A - sqrt(A))

A*B - C*D = (A-K)*(B-K) - (C-K)*(D-K) + K*(A-C+B-D); // Avoid overflow.

Why?

(A-K)*(B-K) = A*B - K*(A+B) + K^2
(C-K)*(D-K) = C*D - K*(C+D) + K^2

=>
(A-K)*(B-K) - (C-K)*(D-K) = A*B - K*(A+B) + K^2 - {C*D - K*(C+D) + K^2}
(A-K)*(B-K) - (C-K)*(D-K) = A*B - C*D - K*(A+B) + K*(C+D) + K^2 - K^2
(A-K)*(B-K) - (C-K)*(D-K) = A*B - C*D - K*(A+B-C-D)

=>
A*B - C*D = (A-K)*(B-K) - (C-K)*(D-K) + K*(A+B-C-D)

=>
A*B - C*D = (A-K)*(B-K) - (C-K)*(D-K) + K*(A-C+B-D)

Note that Because A, B, C and D are big numbers, thus A-C and B-D are small numbers.

참고URL : https://stackoverflow.com/questions/13237046/how-to-avoid-overflow-in-expr-a-b-c-d

반응형