728 x 90

Przekleństwo Kahana: Tanie zabawy floatami

Przekleństwo Kahana: Tanie zabawy floatami

Post dotyczący poprawnego sumowania miliona floatów przypomniał mi przygodę Lorenza z obliczeniami meteo (motyl Lorenza) oraz kilka uwag Nicholasa Highama w książce ,,Accuracy and Stability of Numerical Algorithms” na temat macierzy, które są źle uwarunkowane (ill-conditioned; trochę jak straszliwy karzeł Rumpelstiltskin, ale o nim pod koniec). Marzenia o bezmyślnym stosowaniu floatów pryskają. Na pewno?

Całkiem niedawno natknąłem się na propozycję nowego standardu zapisu liczb, konkurencyjnego wobec floatów, o łacińsko brzmiącej nazwie unums. (Krótki artykuł propagandowy ,,It’s Time for Unums – an Alternative to IEEE 754 Floats and Doubles” autorstwa T. Risse brzmi zachęcająco). Unumsy wymyślił J.L. Gustafson a ich ideę gruntownie przedyskutował w swym książkowym manifeście z 2015 roku o pełnym nadziei tytule ,,The End of Error: Unum Computing”. Nie trzeba było długo czekać na gniew Khana… Znaczy… na ostrą krytykę pomysłu ze strony Kahana – twórcy standardu IEEE 754. Nie wszyscy się jednak wystraszyli, np. twórcy języka Julia chętnie przygarnęli inicjatywę wdrożenia unumsów.

A jakież to przekleństwo Kahana może dotknąć doraźnego kodera? To wyzwania Kahana. Niczym niewierny Tomasz zacząłem testować zapowiadane dziwy w ogólnodostępnym systemie algebry komputerowej Maxima CAS. W zasadzie wszystko się potwierdziło. Z wyjątkiem czterech jedynek Kahana – straszono czterema zerami, a wyszły liczby bliskie 1. Gustafson twierdzi, że jego unumsom testy Kahana niestraszne.

Kody w Maximie testów typu Kahana wylistowano poniżej. (Jeśli ktoś nie ma/nie chce mieć Maximy, to może się też pobawić online w serwisie http://maxima-online.org/. Żeby nie było – zespół od języka Julia również dostarcza zabaw online). Nazwy testów za Gustafsonem.

1. Typowe wyzwanie

E(z):= if (z=0) then 0 else (exp(z)-1)/z$
Q(x):= abs(x-(x^2+1)^(1/2)) - 1/(x+(x^2+1))$
H(x):= E((Q(x))^2)$
'E(z)=E(z); 
'Q(z)=Q(z); 
'H(z)=H(z);
/* winno być: [1;1;1;1] */
/* wychodzi: jak poniżej */
[H(15.0), H(16.0), H(17.0), H(9999.0)];

2. Gładziutka niespodzianka

f(x):= log(abs(3*(1-x)+1))/80 + x^2+1$ /* na przedziale [0.8, 2.0] */
'f(x)=f(x);
/* winno być: -oo w x=4/3 (zwisający jęzor) */
/* wychodzi: minimum w x=0.8 */
wxplot2d([f(x)], [x,0.8,2.0])$

3. Wyzwanie Gustafsona (rzucone Kahanowi)

5.9604644775390625^0.875 = 4.76837158203125; 
/* winno być: to co z prawej */
/* wychodzi: to co z lewej */;

4. Iście królewski ból Kadłubka (czyli młodego Rumpelstiltskina)

P(x,y):= 333.75*y^6+x^2*(11*x^2*y^2-y^6-121*y^4-2) + 5.5*y^8+x/(2*y)$ 
'P(x,y) = P(x,y);
/* winno być: P(x,y) = -0.82739605994682136... (ujemne!) przy x=77617, y=33096 */
/* wychodzi: P(x,y) duża liczba dodatnia */
'P(77617,33096) = P(77617,33096);

Dla leniwych załączam gotowy plik pdf z wynikami: IEEEchallengeKAHAN-Friwebpl2g

A do Rumpelstiltskina może jeszcze kiedyś wrócimy przy okazji algorytmicznej (nieprobabilistycznej) losowości.

Źródło obrazu: opracowanie własne na podstawie zdjęcia z Wikipedii

Leave a Comment

Your email address will not be published. Required fields are marked with *

Cancel reply

Inne artykuły