Względne porównanie liczb zmiennoprzecinkowych

Mam funkcję numeryczną f(x, y)zwracającą podwójną liczbę zmiennoprzecinkową, która implementuje pewną formułę i chcę sprawdzić, czy jest ona poprawna w stosunku do wyrażeń analitycznych dla wszystkich kombinacji parametrów xi yże jestem zainteresowany. Jaki jest właściwy sposób porównania obliczonego i analityczne liczby zmiennoprzecinkowe?

Powiedzmy, że dwie liczby to ai b. Do tej pory upewniłem się, że zarówno błędy bezwzględne ( abs(a-b) < eps), jak i względne ( abs(a-b)/max(abs(a), abs(b)) < eps) są mniejsze niż eps. W ten sposób wykryje niedokładności liczbowe, nawet jeśli liczby będą powiedzmy około 1e-20.

Jednak dzisiaj odkryłem problem, wartość liczbowa ai wartość analityczna bbyły następujące:

In [47]: a                                                                     
Out[47]: 5.9781943146790832e-322

In [48]: b                                                                     
Out[48]: 6.0276008792632078e-322

In [50]: abs(a-b)                                                              
Out[50]: 4.9406564584124654e-324

In [52]: abs(a-b) / max(a, b)                                                  
Out[52]: 0.0081967213114754103

Zatem błąd bezwzględny [50] jest (oczywiście) mały, ale błąd względny [52] jest duży. Pomyślałem więc, że mam błąd w moim programie. Podczas debugowania zdałem sobie sprawę, że liczby te są normalne . Jako taki napisałem następującą procedurę, aby wykonać właściwe porównanie względne:

real(dp) elemental function rel_error(a, b) result(r)
real(dp), intent(in) :: a, b
real(dp) :: m, d
d = abs(a-b)
m = max(abs(a), abs(b))
if (d < tiny(1._dp)) then
    r = 0
else
    r = d / m
end if
end function

Gdzie tiny(1._dp)zwraca 2.22507385850720138E-308 na moim komputerze. Teraz wszystko działa i po prostu dostaję 0 jako błąd względny i wszystko jest w porządku. W szczególności powyższy błąd względny [52] jest błędny, jest po prostu spowodowany niewystarczającą dokładnością liczb denormalnych. Czy moja implementacja rel_errorfunkcji jest poprawna? Czy powinienem po prostu sprawdzić, czy abs(a-b)jest mniejszy niż mały (= denormalny) i zwrócić 0? Czy powinienem sprawdzić inną kombinację, na przykład max(abs(a), abs(b))?

Chciałbym tylko wiedzieć, jaki jest „właściwy” sposób.

Można bezpośrednio sprawdzić denormals korzystających isnormal()z math.h(C99 lub nowszym POSIX.1 lub później). W module Fortran, jeśli moduł ieee_arithmeticjest dostępny, możesz użyć ieee_is_normal(). Aby być bardziej precyzyjnym na temat rozmytej równości, musisz wziąć pod uwagę zmiennoprzecinkową reprezentację denormali i zdecydować, co masz na myśli, aby wyniki były wystarczająco dobre.

Co więcej, aby uwierzyć, że którykolwiek z wyników jest poprawny, musisz mieć pewność, że nie straciłeś zbyt wielu cyfr na etapie pośrednim. Obliczenia z wartościami denormalnymi są generalnie zawodne i należy ich unikać poprzez wewnętrzną zmianę skali algorytmu. Aby upewnić się, że skalowanie wewnętrzne zakończyło się powodzeniem, zalecamy aktywację wyjątków zmiennoprzecinkowych za feenableexcept()pomocą C99 lub ieee_arithmeticmodułu w Fortran.

Mimo że aplikacja może przechwytywać sygnał podniesiony w wyjątkach zmiennoprzecinkowych, wszystkie jądra, które próbowałem zresetować flagę sprzętową, fetestexcept()nie zwracają użytecznego wyniku. Po uruchomieniu z -fp_trapprogramami PETSc (domyślnie) wydrukuje ślad stosu, gdy podniesiony zostanie błąd zmiennoprzecinkowy, ale nie zidentyfikuje linii obrażającej. Jeśli uruchomisz debugger, debuger zachowa flagę sprzętową i zepsuje obrażające wyrażenie. Możesz sprawdzić dokładny powód, dzwoniąc fetestexceptz debuggera, gdzie wynik jest bitowy LUB z następujących flag (wartości mogą się różnić w zależności od komputera, patrz fenv.h; te wartości dotyczą x86-64 z glibc).

• FE_INVALID = 0x1
• FE_DIVBYZERO = 0x4
• FE_OVERFLOW = 0x8
• FE_UNDERFLOW = 0x10
• FE_INEXACT = 0x20

Donald Knuth ma propozycję algorytmu porównywania zmiennoprzecinkowego w tomie 2 „Algorytmy seminumeryczne” z „The Art of Computer Programming”. Został zaimplementowany w C przez Th. Belding (patrz pakiet fcmp ) i jest dostępny w GSL .

Optymalnie zaokrąglone liczby zdormalizowane mogą rzeczywiście mieć wysoki błąd względny. (Spłukanie go do zera, wciąż nazywając go błędem względnym, wprowadza w błąd).

Ale bliskie zeru, obliczanie względnych errosów jest bez znaczenia.

Dlatego nawet przed osiągnięciem liczb zdormalizowanych powinieneś prawdopodobnie przejść do absolutnej dokładności (mianowicie tej, którą chcesz zagwarantować w tym przypadku).

$y$$x$$|y-x|\le absacc+relacc*\max(|x|,|y|)$

Dzięki temu użytkownicy Twojego kodu dokładnie wiedzą, ile naprawdę mają dokładności.