W dniu wtorek, 26 marca 2013 19:02:08 UTC+1 użytkownik Adam Klobukowski napisał:


> Każdy procesor który to implementuje (a ogromna większość implementacji
> matematyki zmiennoprzecinkowej implementuje ten standard) musi liczyć
> dokłądnie tak samo. Jak nie liczy tak samo, to masz np. słynny Pentium bug.


Poniższy program na platformach sprzętowych jakie miałem pod ręką
zadziałał identycznie. Gdy wysłałem znajomemu żeby skompilował i
uruchomił u siebie, wyniki były inne. Ciekawe z czego to wynika.

Implementacja tego drobiazgu w kompilatorze nie wydaje się
trudna, a póki co wszystko wskazuje na to, że winę za niedokładne
obliczenia ponosi kompilator. Jeśli nie jest trudna, to czemu
czasami traci dokładność po małej ilości pętli?


int main( int argc, char *argv[] ) {
const double start = atof(argv[1]);
const double stop = atof(argv[2]);
const double p = atoi(argv[3]);

for( double x=start ; x<=stop ; x+=1.0/1024.0 )
if( pow(x,p) != (x*x*x) ) {
printf("%lf\n",x);
break;
}

return 0;
}

Uruchomienie:
./nazwa_programu 0 1000000 3

Mój wynik:
92681.901367

Szacując na oko, u mnie doszedł daleko, ale innych przerwał dużo
szybciej.

Pozdrawiam

do góry

On Tuesday, 26 March 2013 19:10:04 UTC+1, M.M. wrote:
> W dniu wtorek, 26 marca 2013 19:02:08 UTC+1 użytkownik Adam Klobukowski napisał:
>
> > Nie. W przypadku obliczeń zmiennoprzecinkowych obowiązuje standard IEEE i
> > wszystkie szeroko używane procesory się do niego stosują. Ten standard
> > określa jak procesor ma obliczać i o ile się mylić :)
>
> Była mowa o zwiększaniu dokładności i o ulepszeniach, czyli standard
> zabrania poprawy dokładności?

Standard definiuje miedzy innymi typy danych (np 32, 64, 80 bit float), i w jego
zakresie nie da się zmienić dokładności.

> > Nie. Polecam zapoznać się: http://en.wikipedia.org/wiki/Ieee_floating_point
> > Każdy procesor który to implementuje (a ogromna większość implementacji
> > matematyki zmiennoprzecinkowej implementuje ten standard) musi liczyć
> > dokłądnie tak samo. Jak nie liczy tak samo, to masz np. słynny Pentium bug.
>
> Sprawdzałem pary kompilator/procesor. Właściwie to rzadko zdarzało się, aby
> wyniki były takie same. Jeśli procesory liczą tak samo, to znaczy że
> kompilatory generują różny kod. Będę musiał posprawdzać co się dzieje po
> przeniesieniu wersji binarnej.

Rózne procesory tej samej architektury mogą wykonywać te same operacji tożsamo ale
nie identycznie, choćby dlatego że jedna liczba może mieć kilka reprezentacji
zmiennoprzecinkowych. To się kumuluje w ciągu kilku operacji zmiennoprzecinkowych. Do
niwelowania tego typu problemów (między innymi) służy opcja -march w kompilatorach.
Jak się robi krytyczne obliczenia, to się mierzy błąd operacji zmiennoprzecinkowych i
jak wychodzi za dużo to się kombinuje.

AdamK

do góry

On Tuesday, 26 March 2013 22:00:55 UTC+1, M.M. wrote:
> W dniu wtorek, 26 marca 2013 19:02:08 UTC+1 użytkownik Adam Klobukowski napisał:
>
> > Każdy procesor który to implementuje (a ogromna większość implementacji
> > matematyki zmiennoprzecinkowej implementuje ten standard) musi liczyć
> > dokłądnie tak samo. Jak nie liczy tak samo, to masz np. słynny Pentium bug.
>
> Poniższy program na platformach sprzętowych jakie miałem pod ręką
> zadziałał identycznie. Gdy wysłałem znajomemu żeby skompilował i
> uruchomił u siebie, wyniki były inne. Ciekawe z czego to wynika.
>
> Implementacja tego drobiazgu w kompilatorze nie wydaje się
> trudna, a póki co wszystko wskazuje na to, że winę za niedokładne
> obliczenia ponosi kompilator. Jeśli nie jest trudna, to czemu
> czasami traci dokładność po małej ilości pętli?
>
> int main( int argc, char *argv[] ) {
> const double start = atof(argv[1]);
> const double stop = atof(argv[2]);
> const double p = atoi(argv[3]);
>
> for( double x=start ; x<=stop ; x+=1.0/1024.0 )
> if( pow(x,p) != (x*x*x) ) {
> printf("%lf\n",x);
> break;
> }
>
> return 0;
> }
>
> Uruchomienie:
>
> ./nazwa_programu 0 1000000 3
>
> Mój wynik:
> 92681.901367
>
> Szacując na oko, u mnie doszedł daleko, ale innych przerwał dużo
> szybciej.

U mnie ten sam wynik. Problemem może być choćby 1.0/1024.0 (bo to różnie może
kompilator popracować) no i funkcje atof, atoi że nie wspomnę o rzutowaniu. Oprócz
tego, przy matematyce zmiennoprzecinkowej pow(x,3) wcale nie musi być równe x*x*x.

AdamK

do góry

W dniu wtorek, 26 marca 2013 22:01:28 UTC+1 użytkownik Adam Klobukowski napisał:

> Standard definiuje miedzy innymi typy danych (np 32, 64, 80 bit float), i w
> jego zakresie nie da się zmienić dokładności.
Dla mnie z tego nic nie wynika. Nie wiem co znaczy "zmienić dokładność".

Jeśli typ składa się z iluś tam bitów, jedyną sensowną dokładnością jaką
jestem sobie w stanie wyobrazić, jest dokładność co do najmniej znaczącego
bitu. Jeśli liczba jest sumą całkowitych potęg dwójki i mieści się w zakresie
typu, to nie ma żadnego ważnego powodu aby obliczenia nie były przeprowadzone
dokładnie. Można to osiągnąć i sprzętowo, i software'ow - np. w kompilatorze,
albo w bibliotece do obliczeń zmiennoprzecinkowych. Podobnie sprawa ma
się w każdej innej reprezentacji niż dwójkowa.

Mogę zrozumieć że przyjmuje się różną strategię zaokrągleń ostatniego bitu,
ale jeśli mogą pojawiać się jakieś losowe/deterministyczne i zarazem
niedokładne wartości w najmniej znaczących bitach to sorry, ale dlaczego
mam to nazywać inaczej niż "niedokładne obliczenia tam, gdzie dokładnie
policzyć można"?



> Rózne procesory tej samej architektury mogą wykonywać te same operacji
> tożsamo ale nie identycznie, choćby dlatego że jedna liczba może mieć kilka
> reprezentacji zmiennoprzecinkowych. To się kumuluje w ciągu kilku operacji
> zmiennoprzecinkowych. Do niwelowania tego typu problemów (między innymi)
> służy opcja -march w kompilatorach. Jak się robi krytyczne obliczenia, to się
> mierzy błąd operacji zmiennoprzecinkowych i jak wychodzi za dużo to się
> kombinuje.
Tematyka jest ciekawa, ale byśmy musieli uściślić to o czym rozmawiamy. Ja
trochę się gubię, czasami nie wiem czy w tym wątku rozmawiamy o złym
szkolnictwie, czy o procesorach, czy o kompilatorach. Podobno mam się
czegoś nauczyć, strasznie boli mnie to, że nie wiem czego, tym bardziej
że chcę :)

Że błędy mogą się nakładać - to wiem. Z opcją march kombinowałem,
kilka razy i nie zaobserwowałem żadnych różnic - może trafiłem na
specyficzne kompilatory/procesory/obliczenia - nie wiem. Z opcją
fast-math też kombinowałem - tam gdzie sprawdziłem, też nie było
najmniejszych różnic. Jeśli w jakimś przypadku traciłem dokładność
obliczeń, to kombinowałem - czyli robiłem tak jak mówisz. Przywracałem
stabilność jakimś technikami, nie tylko przy pomocy abs(v1-v2) < epsilon.

W to że 2+2==4.0000000001 nie oznacza niedokładnych obliczeń jak niektórzy
(chyba) sugerowali, nie uwierzę na pewno :)

Więc czego powinienem się douczyć?

Pozdrawiam

do góry

On 2013-03-26, M.M. <m...@g...com> wrote:
> W dniu wtorek, 26 marca 2013 22:01:28 UTC+1 użytkownik Adam Klobukowski napisał:
>
>> Standard definiuje miedzy innymi typy danych (np 32, 64, 80 bit float), i w
>> jego zakresie nie da się zmienić dokładności.
> Dla mnie z tego nic nie wynika. Nie wiem co znaczy "zmienić dokładność".

Pozwól, że zacytuję ciebie sprzed dwóch postów:

>>> Była mowa o zwiększaniu dokładności i o ulepszeniach, czyli standard
>>> zabrania poprawy dokładności?

Zdecyduj się, albo pytasz o zmienianie (zwiększanie) dokładności, albo
nie wiesz, co to ma znaczyć.

> Jeśli typ składa się z iluś tam bitów, jedyną sensowną dokładnością jaką
> jestem sobie w stanie wyobrazić, jest dokładność co do najmniej znaczącego
> bitu.

Błąd reprezentacji to połowa tego.

> Mogę zrozumieć że przyjmuje się różną strategię zaokrągleń ostatniego bitu,
> ale jeśli mogą pojawiać się jakieś losowe/deterministyczne i zarazem
> niedokładne wartości w najmniej znaczących bitach to sorry, ale dlaczego
> mam to nazywać inaczej niż "niedokładne obliczenia tam, gdzie dokładnie
> policzyć można"?

Bo dla takiego brzegowego, rzadko spotykanego przypadku może nie opłacać
się tworzyć wyjątki? Gdy schodzimy do liczb zmiennoprzecinkowych, to
najważniejsza jest możliwość oszacowania błędu obliczeniowego. Dokładne
odwzorowanie liczb z pewnego specyficznego przypadku, który prawie na
pewno nie będzie zachodzić, nie będzie w czołówce problemów do
rozwiązania.

--
Secunia non olet.
Stanislaw Klekot

do góry

Użytkownik "Adam Klobukowski" <a...@g...com> napisał:

> Standard definiuje miedzy innymi typy danych (np 32, 64, 80 bit float),
> i w jego zakresie nie da się zmienić dokładności.

Jaki standard to definiuje ?

AK

do góry

W dniu 2013-03-26 12:25, M.M. pisze:
> Niczego to nie zmienia. Logarytm tez mozna policzyc na mantysie o 2 bity
> dluzszej (np. na precyzji 66bitow) i wynik bylby dokladny.


:)

Szkolny przykład.
x=1/4;
x <- 4*x*(1-x).

Powinien wyjść cykl 1/4, 3/4, 1/4...


for (int i=0;i<1000;i++)
x = 4*x*(1-x);

#include <cstdio>
#include <cmath>
#include <time.h>
#include <limits>

double bla(double x,const int n)
{
for (int i=0;i<n;i++)
{
x= 4*x*(1-x);
}
return x;
}

int main()
{
double a1 = 1.0/4.0;
double a2 = a1+a1*std::numeric_limits<double>::epsilon();
double b1,b2;

for (int n=0; n<=70;n+=10)
{
b1= bla(a1,n); b2= bla(a2,n);
printf("%d %lf %lf %lg\n", n, b1, b2, (b1-b2) );
}
return 0;
}


0 0.250000 0.250000 -5.55112e-017
10 0.750000 0.750000 1.13687e-013
20 0.750000 0.750000 1.16415e-010
30 0.750000 0.750000 1.19209e-007
40 0.750000 0.749878 0.00012208
50 0.750000 0.616384 0.133616
60 0.750000 0.613975 0.136025
70 0.750000 0.995892 -0.245892

Ani dwa, ani nawet 2000 bitów nic nie pomoże.
Najmniejsza odchyłka i trajektorie geometrycznie
się rozbiegają.

Dla lepszego związku z tematem dodam, że dwa bity więcej
nie spowodują, że zawsze dostaniesz dobry wynik. Dwa
bity więcej oznaczają, że najczęściej dostaniesz poprawny
wynik.

Któryś logarytm którejś liczby dobrze zapisanej
w systemie dwójkowym zapisze się źle mimo tych
dodatkowych bitów i katastrofa gotowa.


Zmieńmy naszą iterację na
x= exp(log(4.0)+log(x)+log(1-x));

0 0.250000 0.250000 -5.55112e-017
10 0.750000 0.750000 8.52651e-014
20 0.750000 0.750000 8.73115e-011
30 0.750000 0.750000 8.9407e-008
40 0.749969 0.749878 9.15621e-005
50 0.718126 0.616384 0.101742
60 0.950183 0.613966 0.336216
70 0.840387 0.994657 -0.15427

Odpłynęliśmy równie daleko, mimo, że zarówno argument
naszego wyrażenia, jak i wynik były zapisywalne dokładnie.

pzdr
bartekltg

do góry

> Zdecyduj się, albo pytasz o zmienianie (zwiększanie) dokładności, albo
> nie wiesz, co to ma znaczyć.
Własciwie to już nic nie wiem :) Kilka postów wyżej kazano mi się
czegoś nauczyć, a ja nie mogę, bo nie wiem czego.
Pozdrawiam

do góry

On Tuesday, 26 March 2013 22:53:06 UTC+1, AK wrote:
> U�ytkownik "Adam Klobukowski" <a...@g...com> napisa�:
>
> > Standard definiuje miedzy innymi typy danych (np 32, 64, 80 bit float),
> > i w jego zakresie nie da si� zmieni� dok�adno�ci.
>
> Jaki standard to definiuje ?

http://en.wikipedia.org/wiki/IEEE_floating_point

AdamK

do góry

W dniu wtorek, 26 marca 2013 23:01:23 UTC+1 użytkownik bartekltg napisał:
> W dniu 2013-03-26 12:25, M.M. pisze:
>
> > Niczego to nie zmienia. Logarytm tez mozna policzyc na mantysie o 2 bity
> > dluzszej (np. na precyzji 66bitow) i wynik bylby dokladny.
> :)
> Szkolny przykład.
> x=1/4;
> x <- 4*x*(1-x).
> Powinien wyjść cykl 1/4, 3/4, 1/4...
Mnie chodziło jeszcze o coś innego, ale świetle tego, że procesory
liczą zawsze tak samo, to już jest nieważne :)


> Dla lepszego związku z tematem dodam, że dwa bity więcej
> nie spowodują, że zawsze dostaniesz dobry wynik. Dwa
> bity więcej oznaczają, że najczęściej dostaniesz poprawny
> wynik.
No tak, ale to jednak inny przykład i obliczenia pomimo że
proste, to znacznie bardziej skomplikowane od tych, o jakie
mnie chodziło. Jeśli mamy N liczb i każda z nich jest sumą
całkowitych potęg dwójki, wszystkie są dodatnie, razem
sumują się (dokładnie) do jedynki i mieszczą się w typie, to
w wyniku poniższych operacji na tych liczbach:
i = rand( 1 , N )
j = rand( 1 , N )
tmp = x[i] * 0.5;
x[i] -= tmp;
x[j] += tmp;
Suma nie powinna odbiegać od jeden? Odpalam jeden tego
typu program na dobę i po dobie obliczeń nie mam straty
dokładności.

> Odpłynęliśmy równie daleko, mimo, że zarówno argument
> naszego wyrażenia, jak i wynik były zapisywalne dokładnie.
Hmmm, ale wszystkie pośrednie wyniki też były dokładne? Bo
w przykładzie o jaki mnie chodzi, wszystkie pośrednie wyniki
przynajmniej mogą być dokładne.


Ja idę spać, już naprawdę nic z tego nie rozumiem :)
Pozdrawiam

do góry