slawek
2012-11-01 10:15:13 UTC
Tzw. maszynowy epsilon (see Wikipedia) wynosi nie więcej niż 1.111E-016 dla
liczb 64-bitowych. Taki wynik łatwo otrzymać nawet naiwnym algorytmem, w
którym po kolei sprawdzane są w pętli kolejne wartości epsilon - każda
kolejna nieco (o ułamek procenta) mniejsza od poprzedniej. Algorytm "fast"
adaptacyjnie zmienia krok itd. - nie ma to znacznego wypływu na wynik, ale
liczba kroków jest znacznie mniejsza.
Jednak zaglądając do float.h w MS VS C++ można znaleźć definicję
DBL_EPSILON, wraz ze stosownym komentarzem, 2.22044604925031310000E-016.
Jest to niemal 2 razy więcej, niż naprawdę wynosi epsilon (obliczony właśnie
programem skompilowanym w MSVS C++). "This is not a bug, this is
inaccuracy" - chciałoby się powiedzieć.
Zaglądamy dalej - Matlab - tak ostatnio chwalony - ma wbudowaną funkcję
eps - zgadnijcie co zwraca eps jako wynik liczbowy? Tak, też się zdziwiłem -
przecież Matlab to Matlab.
Jeszcze raz rzut oka do Wikipedii - jest sobie wyraźnie dobra wartość
epsilona dla double w tabelce - ale już np. program w Phytonie i wyniki z
niego - znowu błędne 2.22E-16 . I nie jest to "wina Phytona" - ale po prostu
błąd w programie.
"Phytonowcy", staff MS i ludzie z MathWorks popełnili jeden i ten sam błąd -
dzielili przez dwa. Ciąg wartości x[n], jakie otrzymywali, dla dostatecznie
dużego n nie spełniał nierówności 1.0+x[n] > 1.0. Nie jest źle... jeżeli
pamięta się, że dokładność tak wyznaczonego epsilona wynosi plus minus 50%.
To nawet w większości praktycznych zastosowań wystarcza. Ale nie jest dobrym
pomysłem, by tak niedokładną wartość wrzucać jako wzorcową do float.h - bo
99.8% ludzi będzie w ciemno ufało w nieomylność MS - zwłaszcza, że podane
jest to jako, cyt.:
#define DBL_EPSILON 2.2204460492503131e-016 /* smallest such that
1.0+DBL_EPSILON != 1.0 */
a to sugeruje poprawność wszystkich zapisanych cyfr znaczących. Tymczasem
eps znaleziony przez wykonywanie obliczeń (można oszacować epsilon przez
zapisane 1.0 oraz 1.0+epsilon bit po bicie mantysa i wykładnik - vide
IEEE753) leży gdzieś pomiędzy podanymi zakresami:
naive (no. of steps=36736783):
eps > 1.11022213668763790000E-016
eps <= 1.11022324691088480000E-016
fast (no. of steps=187):
eps > 1.11022302462515650000E-016
eps <= 1.11022302462515680000E-016
slawek
liczb 64-bitowych. Taki wynik łatwo otrzymać nawet naiwnym algorytmem, w
którym po kolei sprawdzane są w pętli kolejne wartości epsilon - każda
kolejna nieco (o ułamek procenta) mniejsza od poprzedniej. Algorytm "fast"
adaptacyjnie zmienia krok itd. - nie ma to znacznego wypływu na wynik, ale
liczba kroków jest znacznie mniejsza.
Jednak zaglądając do float.h w MS VS C++ można znaleźć definicję
DBL_EPSILON, wraz ze stosownym komentarzem, 2.22044604925031310000E-016.
Jest to niemal 2 razy więcej, niż naprawdę wynosi epsilon (obliczony właśnie
programem skompilowanym w MSVS C++). "This is not a bug, this is
inaccuracy" - chciałoby się powiedzieć.
Zaglądamy dalej - Matlab - tak ostatnio chwalony - ma wbudowaną funkcję
eps - zgadnijcie co zwraca eps jako wynik liczbowy? Tak, też się zdziwiłem -
przecież Matlab to Matlab.
Jeszcze raz rzut oka do Wikipedii - jest sobie wyraźnie dobra wartość
epsilona dla double w tabelce - ale już np. program w Phytonie i wyniki z
niego - znowu błędne 2.22E-16 . I nie jest to "wina Phytona" - ale po prostu
błąd w programie.
"Phytonowcy", staff MS i ludzie z MathWorks popełnili jeden i ten sam błąd -
dzielili przez dwa. Ciąg wartości x[n], jakie otrzymywali, dla dostatecznie
dużego n nie spełniał nierówności 1.0+x[n] > 1.0. Nie jest źle... jeżeli
pamięta się, że dokładność tak wyznaczonego epsilona wynosi plus minus 50%.
To nawet w większości praktycznych zastosowań wystarcza. Ale nie jest dobrym
pomysłem, by tak niedokładną wartość wrzucać jako wzorcową do float.h - bo
99.8% ludzi będzie w ciemno ufało w nieomylność MS - zwłaszcza, że podane
jest to jako, cyt.:
#define DBL_EPSILON 2.2204460492503131e-016 /* smallest such that
1.0+DBL_EPSILON != 1.0 */
a to sugeruje poprawność wszystkich zapisanych cyfr znaczących. Tymczasem
eps znaleziony przez wykonywanie obliczeń (można oszacować epsilon przez
zapisane 1.0 oraz 1.0+epsilon bit po bicie mantysa i wykładnik - vide
IEEE753) leży gdzieś pomiędzy podanymi zakresami:
naive (no. of steps=36736783):
eps > 1.11022213668763790000E-016
eps <= 1.11022324691088480000E-016
fast (no. of steps=187):
eps > 1.11022302462515650000E-016
eps <= 1.11022302462515680000E-016
slawek