Jak uniknąć przepełnienia w wyrażeniu. A * B - C * D

Muszę obliczyć wyrażenie, które wygląda następująco:, A*B - C*Dgdzie ich typy to: signed long long int A, B, C, D; Każda liczba może być naprawdę duża (nie przepełniać swojego typu). Chociaż A*Bmoże spowodować przepełnienie, w tym samym czasie wyrażenieA*B - C*D może być naprawdę małe. Jak mogę to poprawnie obliczyć?

Na przykład:, MAX * MAX - (MAX - 1) * (MAX + 1) == 1gdzie in MAX = LLONG_MAX - n- jakaś liczba naturalna.

Wydaje mi się, że to zbyt trywialne. Ale A*Bto ten, który może się przepełnić.

Możesz wykonać następujące czynności bez utraty precyzji

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)

Ten rozkład można przeprowadzić dalej .
Jak zauważył @Gian, może być konieczne zachowanie ostrożności podczas operacji odejmowania, jeśli typ jest długi bez znaku.

Na przykład w przypadku, którego dotyczy pytanie, wystarczy jedna iteracja,

 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

Najprostszym i najbardziej ogólnym rozwiązaniem jest użycie reprezentacji, która nie może się przepełnić, albo za pomocą biblioteki długich liczb całkowitych (np. Http://gmplib.org/ ), albo reprezentując za pomocą struktury lub tablicy i implementując rodzaj długiego mnożenia ( tj. rozdzielenie każdej liczby na dwie 32-bitowe połowy i wykonanie mnożenia jak poniżej:

(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

Zakładając, że wynik końcowy mieści się w 64 bitach, tak naprawdę nie potrzebujesz większości bitów R3 i żadnego R4

Zauważ, że nie jest to standardowe, ponieważ opiera się na zawijanym przepełnieniu ze znakiem. (GCC ma flagi kompilatora, które to umożliwiają.)

Ale jeśli wykonasz wszystkie obliczenia w programie long long, wynik bezpośredniego zastosowania wzoru:
(A * B - C * D)będzie dokładny, o ile poprawny wynik będzie pasował do a long long.

Oto obejście, które polega tylko na zachowaniu zdefiniowanym przez implementację polegającym na rzutowaniu liczby całkowitej bez znaku na liczbę całkowitą ze znakiem. Ale dziś można oczekiwać, że zadziała to na prawie każdym systemie.

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

To rzutuje dane wejściowe do unsigned long longmiejsca, w którym zachowanie przepełnienia jest gwarantowane jako zawijane przez standard. Rzutowanie z powrotem na podpisaną liczbę całkowitą na końcu jest częścią zdefiniowaną przez implementację, ale dziś będzie działać w prawie wszystkich środowiskach.

Jeśli potrzebujesz bardziej pedantycznego rozwiązania, myślę, że musisz użyć „długiej arytmetyki”

To powinno działać (myślę):

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

Oto moje wyprowadzenie:

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)

Możesz rozważyć obliczenie największego wspólnego współczynnika dla wszystkich wartości, a następnie podzielenie ich przez ten współczynnik przed wykonaniem operacji arytmetycznych, a następnie ponowne pomnożenie. Zakłada się, że taki czynnik istnieje jednak (na przykład w przypadku A, B, CiD stało się względnie pierwsze, nie będą mieć wspólny czynnik).

Podobnie możesz rozważyć pracę na skalach logarytmicznych, ale będzie to trochę przerażające, z zastrzeżeniem dokładności numerycznej.

Jeśli wynik mieści się w długim długim int, to wyrażenie A * BC * D jest w porządku, ponieważ wykonuje arytmetyczną mod 2 ^ 64 i da prawidłowy wynik. Problem polega na tym, aby wiedzieć, czy wynik mieści się w długiej, długiej int. Aby to wykryć, możesz użyć następującej sztuczki, używając podwójnych:

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

Problem z tym podejściem polega na tym, że jesteś ograniczony precyzją mantysy podwójnych (54 bity?), Więc musisz ograniczyć iloczyn A * B i C * D do 63 + 54 bitów (lub prawdopodobnie trochę mniej).

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

następnie

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

Możesz zapisać każdą liczbę w tablicy, przy czym każdy element jest cyfrą i wykonywać obliczenia jako wielomiany . Weź wynikowy wielomian, który jest tablicą, i oblicz wynik, mnożąc każdy element tablicy przez 10 do potęgi pozycji w tablicy (pierwsza pozycja jest największa, a ostatnia zero).

Liczbę 123można wyrazić jako:

123 = 100 * 1 + 10 * 2 + 3

dla którego po prostu tworzysz tablicę [1 2 3].

Robisz to dla wszystkich liczb A, B, C i D, a następnie mnożysz je jako wielomiany. Gdy masz wynikowy wielomian, po prostu rekonstruujesz z niego liczbę.

Podczas gdy jeden signed long long intnie wytrzyma A*B, dwóch z nich to zrobi. A*BMożna więc rozłożyć na trzy wyrażenia o różnych wykładnikach, z których każdy pasuje do jednego 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;

To samo dotyczy C*D .

Idąc prostą drogą, można było dokonać subrakcji na każdej parze AB_ii CD_ipodobnie, używając dodatkowego bitu przenoszenia (dokładnie 1-bitowej liczby całkowitej) dla każdej. Więc jeśli powiemy E = A * BC * D, otrzymasz coś takiego:

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

Kontynuujemy, przenosząc górną połowę E_10do E_20(przesuń o 32 i dodaj, a następnie usuń górną połowę E_10).

Teraz możesz pozbyć się bitu przenoszonego E_11, dodając go z odpowiednim znakiem (uzyskanym z części nienoszonej) do E_20. Jeśli spowoduje to przepełnienie, wynik również nie będzie pasował.

E_10teraz ma wystarczająco dużo miejsca, aby wziąć górną połowę z E_00 (przesuń, dodaj, wymaż) i przenoszoną część E_01.

E_10może być teraz ponownie większy, więc powtarzamy transfer do E_20.

W tym momencie E_20musi wynosić zero, w przeciwnym razie wynik nie będzie pasował. W górnej połowieE_10 jest również pusta w wyniku przeniesienia.

Ostatnim krokiem jest ponowne przeniesienie dolnej połowy E_20do E_10.

Jeśli oczekiwanie, E=A*B+C*Dktóre pasowałoby do signed long long intblokad, mamy teraz

E_20=0
E_10=0
E_00=E

Jeśli wiesz, że wynik końcowy można przedstawić w typie liczby całkowitej, możesz szybko wykonać to obliczenie, korzystając z poniższego kodu. Ponieważ standard C określa, że ​​arytmetyka bez znaku jest arytmetyką modulo i nie powoduje przepełnienia, można użyć typu bez znaku do wykonania obliczeń.

Poniższy kod zakłada, że ​​istnieje typ bez znaku o tej samej szerokości i że typ ze znakiem wykorzystuje wszystkie wzorce bitowe do reprezentowania wartości (bez reprezentacji pułapek, minimum typu ze znakiem jest ujemną wartością połowy modułu typu bez znaku). Jeśli tak się nie stanie w implementacji C, można w tym celu wprowadzić proste poprawki w procedurze ConvertToSigned.

Poniższe zastosowania signed chari unsigned charzademonstrowanie kodu. Na potrzeby swojej implementacji zmień definicję Signedto typedef signed long long int Signed;i definicję Unsignedto 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;
}

Możesz spróbować rozbić równanie na mniejsze komponenty, które się nie przepełniają.

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

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

    where K = B - N
          J = D - M

Jeśli komponenty nadal się przepełniają, możesz je rekurencyjnie podzielić na mniejsze komponenty, a następnie ponownie połączyć.

Być może nie omówiłem wszystkich przypadków skrajnych, ani nie przetestowałem tego rygorystycznie, ale implementuje to technikę, którą pamiętam z lat 80., kiedy próbowałem wykonywać 32-bitowe obliczenia liczb całkowitych na 16-bitowym procesorze. Zasadniczo dzielisz 32 bity na dwie 16-bitowe jednostki i pracujesz z nimi osobno.

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

  }

}

Wydruki:

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

co wydaje mi się, że działa.

Założę się, że przegapiłem niektóre subtelności, takie jak obserwowanie przepełnienia znaków itp., Ale myślę, że istota jest tam.

Ze względu na kompletność, ponieważ nikt o tym nie wspomniał, niektóre kompilatory (np. GCC) obecnie dostarczają 128-bitową liczbę całkowitą.

Zatem prostym rozwiązaniem mogłoby być:

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

AB-CD = (AB-CD) * AC / AC = (B/C-D/A)*A*C. Ani, B/Cani nie D/Amogą się przepełniać, więc (B/C-D/A)najpierw oblicz . Ponieważ ostateczny wynik nie przepełni się zgodnie z twoją definicją, możesz bezpiecznie wykonać pozostałe mnożenia i obliczyć, (B/C-D/A)*A*Cktóry wynik jest wymagany.

Uwaga, jeśli twój wkład może być również bardzo mały , B/Club D/Amoże się przepełnić. Jeśli to możliwe, mogą być wymagane bardziej złożone manipulacje zgodnie z kontrolą wejściową.

Wybierz K = a big number(np. K = A - sqrt(A))

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

Czemu?

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

Zauważ, że ponieważ A, B, C i D są dużymi liczbami, A-Ca więc i B-Dsą małymi liczbami.