Solidne obliczenie średniej z dwóch liczb zmiennoprzecinkowych?

Pozwól x, ybyć dwóch liczb zmiennoprzecinkowych. Jaki jest właściwy sposób na obliczenie ich średniej?

Naiwny sposób (x+y)/2może doprowadzić do przepełnienia kiedy xi ysą zbyt duże. Myślę, że 0.5 * x + 0.5 * ymoże lepiej, ale wiąże się to z dwoma multiplikacjami (co może być nieefektywne) i nie jestem pewien, czy jest wystarczająco dobry. Czy jest lepszy sposób?

Innym pomysłem, z którym bawiłem się, jest to, (y/2)(1 + x/y)czy x<=y. Ale znowu nie jestem pewien, jak to przeanalizować i udowodnić, że spełnia moje wymagania.

Ponadto potrzebuję gwarancji, że obliczoną średnią będzie >= min(x,y)i <= max(x,y). Jak wskazano w odpowiedzi Dona Hatcha , może lepszym sposobem postawienia tego pytania jest: Jaka jest implementacja średniej z dwóch liczb, która zawsze daje najbardziej możliwy dokładny wynik? To znaczy, jeśli xi yczy liczby zmiennoprzecinkowe są najbliżej, jak obliczyć liczbę zmiennoprzecinkową (x+y)/2? W takim przypadku obliczona średnia jest automatycznie >= min(x,y)i <= max(x,y). Szczegółowe informacje można znaleźć w odpowiedzi Dona Hatcha .

Uwaga: Moim priorytetem jest solidna dokładność. Wydajność jest zbywalna. Jeśli jednak istnieje wiele niezawodnych i dokładnych algorytmów, wybrałbym najbardziej wydajny.

Myślę, że dokładność i stabilność algorytmów numerycznych Highama dotyczy tego, w jaki sposób można analizować tego rodzaju problemy. Patrz rozdział 2, zwłaszcza ćwiczenie 2.8.

W tej odpowiedzi chciałbym wskazać coś, co tak naprawdę nie zostało poruszone w książce Highama (wydaje się, że nie jest to zbyt powszechnie znane). Jeśli jesteś zainteresowany udowodnieniem właściwości prostych algorytmów numerycznych takich jak te, możesz skorzystać z mocy nowoczesnych solverów SMT ( Teorie satysfakcji modulo ), takich jak z3 , używając pakietu takiego jak sbv w Haskell. Jest to nieco łatwiejsze niż używanie ołówka i papieru.

Załóżmy, że podano mi i chciałbym wiedzieć, czy spełnia . Poniższy kod Haskellz = ( x + y ) / 2 x ≤ z ≤ $0\leq x\leq y$$z=(x+y)/2$$x\leq z\leq y$

import Data.SBV

test1 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test1 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ 0 .<= x &&& x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

test2 :: (SFloat -> SFloat -> SFloat) -> Symbolic SBool
test2 fun =
  do [x, y] <- sFloats ["x", "y"]
     constrain $ bnot (isInfiniteFP x) &&& bnot (isInfiniteFP y)
     constrain $ x .<= y
     let z = fun x y
     return $ x .<= z &&& z .<= y

pozwoli mi to zrobić automatycznie . Tutaj test1 funjest założenie , że dla wszystkich skończonej pływaki z .x , y 0 ≤ x ≤ $x \leq \mathit{fun}(x,y) \leq y$$x,y$$0\leq x\leq y$

λ> prove $ test1 (\x y -> (x + y) / 2)
Falsifiable. Counter-example:
  x = 2.3089316e36 :: Float
  y = 3.379786e38 :: Float

To się przelewa. Załóżmy, że teraz biorę inną formułę:$z=x/2+y/2$

λ> prove $ test1 (\x y -> x/2 + y/2)
Falsifiable. Counter-example:
  x = 2.3509886e-38 :: Float
  y = 2.3509886e-38 :: Float

Nie działa (ze względu na stopniowe niedopełnienie: , co może być nieintuicyjne, ponieważ cała arytmetyka ma wartość base-2).$(x/2)\times2 \neq x$

Teraz spróbuj :$z=x + (y-x)/2$

λ> prove $ test1 (\x y -> x + (y-x)/2)
Q.E.D.

Pracuje! Q.E.D.To dowód , że test1własność zachodzi dla wszystkich pływaków, jak zdefiniowano powyżej.

Co z tym samym, ale ograniczonym do (zamiast )?$x\leq y$$0\leq x\leq y$

λ> prove $ test2 (\x y -> x + (y-x)/2)
Falsifiable. Counter-example:
  x = -3.1300826e34 :: Float
  y = 3.402721e38 :: Float

Okej, więc jeśli przepełni, co powiesz na ?$y-x$$z = x + (y/2-x/2)$

λ> prove $ test2 (\x y -> x + (y/2 - x/2))
Q.E.D.

Wygląda więc na to, że spośród wzorów, które tu wypróbowałem, wydaje się działać (również z dowodem). Metoda solvera SMT wydaje mi się o wiele szybszym sposobem odpowiedzi na podejrzenia dotyczące prostych wzorów zmiennoprzecinkowych niż analizowanie błędów zmiennoprzecinkowych ołówkiem i papierem.$x + (y/2 - x/2)$

Wreszcie cel dokładności i stabilności często stoi w sprzeczności z celem wydajności. Jeśli chodzi o wydajność, tak naprawdę nie widzę, jak możesz sobie radzić lepiej niż , zwłaszcza, że ​​kompilator nadal będzie cię ciężko tłumaczyć, tłumacząc to na instrukcje maszynowe.$(x+y)/2$

PS Wszystko to z arytmetyką zmiennoprzecinkową IEEE754 o pojedynczej precyzji. Sprawdziłem z podwójnej precyzji arytmetyki (wymienić z ) i działa zbyt.$x \leq x + (y/2-x/2) \leq y$SFloatSDouble

PPS Jedną z rzeczy, o których należy pamiętać przy implementacji tego w kodzie jest to, że flagi kompilatora takie jak -ffast-math(niektóre formy takich flag są czasami domyślnie włączone w niektórych popularnych kompilatorach) nie spowodują arytmetyki IEEE754, co unieważni powyższe dowody. Jeśli używasz flag, które umożliwiają np. Optymalizacje dodawania skojarzonego, nie ma sensu robić niczego innego niż .$(x+y)/2$

PPPS Dałam się trochę ponieść spojrzeniu tylko na proste wyrażenia algebraiczne bez warunków warunkowych. Don Hatch „s formuła jest ściśle lepiej.

Po pierwsze, zauważ, że jeśli masz metodę, która daje najdokładniejszą odpowiedź we wszystkich przypadkach, spełni ona wymagany warunek. (Należy pamiętać, że mówię najdokładniejszą odpowiedź zamiast z najdokładniejszej odpowiedzi, ponieważ nie może być dwóch zwycięzców.) Dowód: Jeśli, przeciwnie, trzeba dokładnej-as-możliwa odpowiedź, która ma nie spełniają wymaganego warunku, że oznacza albo (w którym przypadku jest lepsza odpowiedź, sprzeczność), albo (w którym przypadku jest lepsza odpowiedź, sprzeczność).answer<min(x,y)<=max(x,y)min(x,y)min(x,y)<=max(x,y)<answermax(x,y)

Myślę więc, że to oznacza, że ​​twoje pytanie sprowadza się do znalezienia najdokładniejszej możliwej odpowiedzi. Zakładając, że arytmetyka IEEE754 jest w toku, proponuję:

if max(abs(x),abs(y)) >= 1.:
    return x/2. + y/2.
else:
    return (x+y)/2.

Mój argument, że daje to najdokładniejszą odpowiedź, jest dość żmudną analizą przypadku. Tutaj idzie:

• Sprawa max(abs(x),abs(y)) >= 1.:

  • Podklucz ani x, ani y nie są denormalizowane: w tym przypadku obliczona odpowiedź x/2.+y/2.manipuluje tymi samymi mantysami, a zatem daje dokładnie taką samą odpowiedź, jak obliczenie (x+y)/2dałoby, gdybyśmy przyjęli rozszerzone wykładniki, aby zapobiec przepełnieniu. Ta odpowiedź może zależeć od trybu zaokrąglania, ale w każdym przypadku IEEE754 gwarantuje, że jest to najlepsza możliwa odpowiedź (z faktu, że obliczone x+ygwarantuje najlepsze przybliżenie do matematyki x + y, a podział przez 2 jest w tym dokładny walizka).

  • Podtekst x jest zdenormalizowany (i tak abs(y)>=1):

    answer = x/2. + y/2. = y/2. since abs(x/2.) is so tiny compared to abs(y/2.) = the exact mathematical value of y/2 = a best possible answer.

  • Podsekcja y jest zdenormalizowana (i tak abs(x)>=1): analogicznie.

• Sprawa max(abs(x),abs(y)) < 1.:
  • Podliczona obliczona x+yjest albo niez denormalizowana, albo zdenormalizowana, a „- nawet”: Chociaż obliczona x+ymoże nie być dokładna, IEEE754 gwarantuje, że jest najlepszym możliwym przybliżeniem do matematyki x + y. W tym przypadku kolejny podział przez 2 w wyrażeniu (x+y)/2.jest dokładny, więc obliczona odpowiedź (x+y)/2.jest najlepszym możliwym przybliżeniem do matematyki (x + y) / 2.
  • Podklucz obliczonego x+yjest zdenormalizowany i „nieparzysty”: W tym przypadku dokładnie jeden z x, y również musi być zdenormalizowany - i - „nieparzysty”, co oznacza, że ​​drugi z x, y jest zdenormalizowany znakiem przeciwnym, a zatem obliczone x+yjest dokładnie matematyczne x + y, a zatem obliczone (x+y)/2.jest gwarantowane przez IEEE754, aby być najlepszym możliwym przybliżeniem do matematyki (x + y) / 2.

W przypadku binarnych formatów zmiennoprzecinkowych IEEE-754, których przykładem jest binary64obliczenie (podwójnej precyzji), S. Boldo formalnie udowodnił, że przedstawiony poniżej prosty algorytm zapewnia poprawnie zaokrągloną średnią.

Sylvie Boldo, „Formalna weryfikacja programów obliczających średnią zmiennoprzecinkową”. Na międzynarodowej konferencji na temat formalnych metod inżynieryjnych , s. 17–32. Springer, Cham, 2015. ( projekt online )

$(x+y)/2$$x/2 + y/2$binary64$C \in [2^{-967}, 2^{970}]$$C$ tak aby zapewnić najlepszą wydajność dla konkretnego przypadku użycia.

Daje to następujący przykładowy ISO-C99kod:

double average (double x, double y) 
{
    const double C = 1; /* 0x1p-967 <= C <= 0x1p970 */
    return (C <= fabs (x)) ? (x / 2 + y / 2) : ((x + y) / 2);
}

W ostatnich pracach uzupełniających S. Boldo i współautorzy pokazali, jak osiągnąć najlepsze możliwe wyniki dla dziesiętnych formatów zmiennoprzecinkowych IEEE-754, wykorzystując operacje fuzji wielokrotnego dodawania (FMA) i dobrze znaną precyzję podwajanie bloku konstrukcyjnego (TwoSum):

Sylvie Boldo, Florian Faissole i Vincent Tourneur, „Formalnie ustalony algorytm obliczania poprawnej średniej liczb dziesiętnych zmiennoprzecinkowych”. W 25. sympozjum IEEE na temat arytmetyki komputerowej (ARITH 25) , czerwiec 2018 r., S. 69–75. ( projekt online )

Chociaż nie może być super-wydajny wydajność mądry, istnieje bardzo prosty sposób (1) upewnij się, że żaden z tych liczb jest większa niż którakolwiek xlub y(bez przepełnienia) i (2) utrzymać zmiennoprzecinkowych jako „dokładne”, jak możliwe (i (3) , jako dodatkowy bonus, nawet jeśli stosowane jest odejmowanie, żadne wartości nie będą nigdy przechowywane jako liczby ujemne.

float difference = max(x, y) - min(x, y);
return min(x, y) + (difference / 2.0);

W rzeczywistości, jeśli naprawdę chcesz uzyskać dokładność, nie musisz nawet dokonywać podziału na miejscu; po prostu zwróć wartości min(x, y)i differencektórych możesz użyć, aby uprościć logicznie lub manipulować później.

Konwertuj na wyższą precyzję, dodaj tam wartości i przekonwertuj z powrotem.

Wyższa precyzja nie powinna powodować przepełnienia, a jeśli oba są w prawidłowym zakresie zmiennoprzecinkowym, obliczona liczba również powinna być w środku.

I powinna znajdować się pomiędzy nimi, w najgorszym przypadku tylko połowa większej liczby, jeśli prewencja nie jest wystarczająca.

Teoretycznie x/2można go obliczyć, odejmując 1 od mantysy.

Jednak faktyczne implementowanie takich operacji bitowych niekoniecznie jest proste, szczególnie jeśli nie znasz formatu liczb zmiennoprzecinkowych.

Jeśli możesz to zrobić, cała operacja zostaje zredukowana do 3 dodawania / odejmowania, co powinno być znaczącą poprawą.

Myślałem w tym samym kierunku, co @Roland Heath, ale nie mogę jeszcze komentować, oto moje zdanie:

x/2można obliczyć, odejmując 1 od wykładnika (nie mantysy, odejmując 1 od mantysy odejmując 2^(value_of_exponent-length_of_mantissa)od wartości ogólnej).

Załóżmy, że bez ograniczenia ogólnego przypadku x < y. (Jeśli x > y, ponownie oznacz zmienne. Jeśli x = y, (x+y) / 2jest banalne.)

• Przekształć (x+y) / 2w x/2 + y/2, który można wykonać przez dwie odejmowanie liczb całkowitych (jedną od wykładnika)
  • Jednak wykładnik ma dolną granicę w zależności od reprezentacji. Jeśli wykładnik jest już minimalny przed odjęciem 1, ta metoda będzie wymagać specjalnej obsługi przypadków. Minimalny wykładnik włączony xsprawi, że będzie on x/2mniejszy niż reprezentowalny (zakładając, że mantysa jest reprezentowana z ukrytym prowadzeniem 1).
  • Zamiast odejmować 1 od wykładnika x, xprzesuń mantysę w prawo o jeden (i dodaj ukrytą wiodącą 1, jeśli istnieje).
  • Odejmij 1 od wykładnika y, jeśli nie jest minimalny. Jeśli jest minimalna (y jest większe niż x, ze względu na mantysę), przesuń mantysę w prawo o jeden (dodaj ukryte prowadzenie 1, jeśli istnieje).
  • Przesuń nową mantysę xw prawo zgodnie z wykładnikiem y.
  • Wykonaj dodawanie liczb całkowitych na mantissae, chyba że mantysa xzostała całkowicie przesunięta. Jeśli oba wykładniki były minimalne, wiodące przepełniłyby się, co jest w porządku, ponieważ to przepełnienie powinno stać się ponownie domyślnym wiodącym.
• i jeden dodatek zmiennoprzecinkowy.
  • Nie mogę tutaj wymyślić żadnego specjalnego przypadku; z wyjątkiem zaokrąglania, które ma również zastosowanie do przesunięć opisanych powyżej.