Właściwości 80-bitowych rozszerzonych obliczeń precyzyjnych, począwszy od argumentów podwójnej precyzji

Question

Oto dwie implementacje funkcji interpolacji. Argument u1 jest zawsze między 0. i 1..

#include <stdio.h> 

double interpol_64(double u1, double u2, double u3) 
{ 
    return u2 * (1.0 - u1) + u1 * u3; 
} 

double interpol_80(double u1, double u2, double u3) 
{ 
    return u2 * (1.0 - (long double)u1) + u1 * (long double)u3; 
} 

int main() 
{ 
    double y64,y80,u1,u2,u3; 
    u1 = 0.025; 
    u2 = 0.195; 
    u3 = 0.195; 
    y64 = interpol_64(u1, u2, u3); 
    y80 = interpol_80(u1, u2, u3); 
    printf("u2: %a\ny64:%a\ny80:%a\n", u2, y64, y80); 
} 

na ścisłej IEEE 754 platformie z 80-bitowymi long double s, wszystkie obliczenia wykonywane są w interpol_64() według IEEE 754 podwójnej precyzji, aw interpol_80() w 80-bitowej precyzji rozszerzonego. drukuje programowe:

u2: 0x1.8f5c28f5c28f6p-3 
y64:0x1.8f5c28f5c28f5p-3 
y80:0x1.8f5c28f5c28f6p-3 

Jestem zainteresowany w nieruchomości „wynik zwrócony przez funkcję jest zawsze pomiędzy u2 i u3”. Ta właściwość ma wartość false wynoszącą interpol_64(), co pokazują wartości w powyższym zbiorze wartości .

Czy nieruchomość ma szansę być prawdziwa z interpol_80()? Jeśli nie, to co jest kontrprzykładem? Czy to pomaga, jeśli wiemy, że istnieje między nimi minimalna odległość? Czy istnieje metoda określania znaczenia i szerokości dla obliczeń pośrednich, przy których właściwość będzie gwarantowana?

EDYTOWANIE: dla wszystkich wartości losowych, które wypróbowałem, właściwość zatrzymana, gdy obliczenia pośrednie zostały wykonane wewnętrznie w rozszerzonej precyzji. Jeśli interpol_80() zajął argumenty, byłoby stosunkowo łatwo zbudować kontrprzykład, ale tutaj jest pytanie o funkcję, która pobiera argumenty double. To sprawia, że ​​znacznie trudniej jest zbudować kontrprzykład, jeśli taki istnieje.

---

Uwaga: generujące kompilator instrukcje x87 może generować ten sam kod interpol_64() i interpol_80(), ale to jest styczny do mojego pytania.

Answer 1

Tak, interpol_80() jest bezpieczny, pokażmy to.

Problem stanowi, że wejścia są 64Bits unosić

rnd64(ui) = ui 

Wynik jest dokładnie (przy założeniu, * i + są operacje matematyczne)

r = u2*(1-u1)+(u1*u3) 

Optymalna wartość zwracana w zaokrągleniu do 64 bit float jest

r64 = rnd64(r) 

Ponieważ posiadamy te właściwości

u2 <= r <= u3 

gwarantuje, że

rnd64(u2) <= rnd64(r) <= rnd64(u3) 
u2 <= r64 <= u3 

Konwersja do 80bits U1, U2, U3 są zbyt dokładne.

rnd80(ui)=ui 

Teraz załóżmy 0 <= u2 <= u3, następnie występował z niedokładne operacji zmiennoprzecinkowych prowadzi do co najwyżej 4 błędy zaokrągleń:

rf = rnd(rnd(u2*rnd(1-u1)) + rnd(u1*u3)) 

Zakładając zaokrąglenie do najbliższej nawet, to będzie co najwyżej 2 ULP off dokładny wartość. przypadku zaokrąglenia jest wykonywana z 64 bitów lub 80 bitów FLOAT unosi:

r - 2 ulp64(r) <= rf64 <= r + 2 ulp64(r) 
r - 2 ulp80(r) <= rf80 <= r + 2 ulp80(r) 

rf64 może być wyłączone przez 2 ULP tak interpol-64() jest niebezpieczna, ale co rnd64(rf80)?
Możemy powiedzieć, że:

rnd64(r - 2 ulp80(r)) <= rnd64(rf80) <= rnd64(r + 2 ulp80(r)) 

Od 0 <= u2 <= u3, następnie

ulp80(u2) <= ulp80(r) <= ulp80(r3) 
rnd64(u2 - 2 ulp80(u2)) <= rnd64(r - 2 ulp80(r)) <= rnd64(rf80) 
rnd64(u3 + 2 ulp80(u3)) >= rnd64(r + 2 ulp80(r)) >= rnd64(rf80) 

Na szczęście, jak każdy numer w zakresie (u2-ulp64(u2)/2 , u2+ulp64(u2)/2) otrzymujemy

rnd64(u2 - 2 ulp80(u2)) = u2 
rnd64(u3 + 2 ulp80(u3)) = u3 

od ulp80(x)=ulp62(x)/2^(64-53)

Mamy więc uzyskać dowód

u2 <= rnd64(rf80) <= u3 

Dla u2 < = U3 < = 0, możemy zastosować tę samą dowód łatwo.

Ostatni przypadek do zbadania to u2 < = 0 < = u3. Jeśli odejmiemy 2 duże wartości, wynik może wynosić do ulp (duży)/2, a nie ulp (duży-duży)/2 ...
Zatem twierdzenie to zrobiliśmy nie posiada już:

r - 2 ulp64(r) <= rf64 <= r + 2 ulp64(r) 

szczęście u2 <= u2*(1-u1) <= 0 <= u1*u3 <= u3 i ten jest zachowany po zaokrągleniu

u2 <= rnd(u2*rnd(1-u1)) <= 0 <= rnd(u1*u3) <= u3 

Tak więc od ilości dodanych mają przeciwny znak:

u2 <= rnd(u2*rnd(1-u1)) + rnd(u1*u3) <= u3 

to samo, co po zaokrągleniu, więc możemy jeszcze raz zagwarantować

u2 <= rnd64(rf80) <= u3 

QED

Aby być kompletne powinniśmy dbać Brak reprezentacji wejść (stopniowy niedomiar), ale mam nadzieję, że nie będzie to błędne z testów warunków skrajnych. I nie pokaże, co się dzieje z tymi ...

EDIT:

Oto obserwacji jako następujące twierdzenie było nieco przybliżony i generowane kilka uwag, gdy 0 = u2 < < = u3

r - 2 ulp80(r) <= rf80 <= r + 2 ulp80(r) 

Możemy napisać następujące nierówności:

rnd(1-u1) <= 1 
rnd(1-u1) <= 1-u1+ulp(1)/4 
u2*rnd(1-u1) <= u2 <= r 
u2*rnd(1-u1) <= u2*(1-u1)+u2*ulp(1)/4 
u2*ulp(1) < 2*ulp(u2) <= 2*ulp(r) 
u2*rnd(1-u1) < u2*(1-u1)+ulp(r)/2 

Dla następna operacja zaokrąglania, używamy

ulp(u2*rnd(1-u1)) <= ulp(r) 
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)/2 + ulp(u2*rnd(1-u1))/2 
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r)/2 + ulp(r)/2 
rnd(u2*rnd(1-u1)) < u2*(1-u1)+ulp(r) 

Dla drugiej części sumy, mamy:

u1*u3 <= r 
rnd(u1*u3) <= u1*u3 + ulp(u1*u3)/2 
rnd(u1*u3) <= u1*u3 + ulp(r)/2 

rnd(u2*rnd(1-u1))+rnd(u1*u3) < u2*(1-u1)+u1*u3 + 3*ulp(r)/2 
rnd(rnd(u2*rnd(1-u1))+rnd(u1*u3)) < r + 3*ulp(r)/2 + ulp(r+3*ulp(r)/2)/2 
ulp(r+3*ulp(r)/2) <= 2*ulp(r) 
rnd(rnd(u2*rnd(1-u1))+rnd(u1*u3)) < r + 5*ulp(r)/2 

Nie udowodnić pierwotnego roszczenia, ale nie tak daleko ...

Answer 2

Głównym źródłem utraty dokładności w interpol_64 jest multiplikacja. Pomnożenie dwóch 53-bitowych mantyli daje 105- lub 106-bitową (w zależności od tego, czy przenosi się na wysoki bit) mantysę. Jest on zbyt duży, aby zmieścić się w 80-bitowej rozszerzonej wartości dokładności, więc generalnie, w wersji 80-bitowej będzie również występować utrata precyzji. Dokładne oszacowanie, kiedy to się dzieje, jest bardzo trudne; Najłatwiej powiedzieć, że dzieje się tak, gdy gromadzą się błędy zaokrąglania. Zwróć uwagę, że przy dodawaniu dwóch terminów jest również mały krok zaokrąglania.

Większość ludzi pewnie tylko rozwiązać ten problem za pomocą funkcji takich jak:

double interpol_64(double u1, double u2, double u3) 
{ 
    return u2 + u1 * (u3 - u2); 
} 

Ale wygląda na to, czego szukasz wgląd kwestii zaokrąglania, a nie lepszego wdrożenia.