Mamy dane wejściowe:
tab[N] - tablica liczb zmiennoprzecinkowych

Dane wyjsciowe:
Iloczyn tab[0] * tab[1] * ... * tab[N-1]

Pytanie:
Jak posortować dane w tab aby mnożenie było maksymalnie stabilne
numerycznie?

Przy dodawaniu najlepiej posortować od najmniejszych co do wartości
bezwzględnej. A przy mnożeniu ?

Pozdrawiam

do góry

On 12.10.2016 17:41, M.M. wrote:
>
> Mamy dane wejściowe:
> tab[N] - tablica liczb zmiennoprzecinkowych
>
> Dane wyjsciowe:
> Iloczyn tab[0] * tab[1] * ... * tab[N-1]
>
> Pytanie:
> Jak posortować dane w tab aby mnożenie było maksymalnie stabilne
> numerycznie?

Nie ma znaczenia.
Mnożysz mantysy, które zawsze są w przedziale [0.5,1)
cechy dodajesz stałoprzecinkowo.
["Ty w sensie komputer jak mnożysz zmienne float/double",
nie trzeba nic ręcznie poprawiać].


> Przy dodawaniu najlepiej posortować od najmniejszych co do wartości
> bezwzględnej.

Nie jest to prawda.

Jest to jedna z metod sumowania, ale okazuje się, że wcale
nie najlepsza (najlepsza w sensie da najlepszy wynik
dla każdego ciągu).
Gdzieś kiedyś miałęm dobrą pracę na ten temat,
a do tego napisany zestaw kodów dowolnej pracyzji
do takiego sumowania (od najmniejszej od największej,
ze stosem, Kahanem, parami...)
ale za szybko tego nie odkopię.

pzdr
bartekltg

do góry

On Wednesday, October 12, 2016 at 6:42:55 PM UTC+2, bartekltg wrote:
> Nie ma znaczenia.
> Mnożysz mantysy, które zawsze są w przedziale [0.5,1)
> cechy dodajesz stałoprzecinkowo.
> ["Ty w sensie komputer jak mnożysz zmienne float/double",
> nie trzeba nic ręcznie poprawiać].

Fajny sposób.
Pozdrawiam

do góry

On 12.10.2016 19:55, M.M. wrote:
> On Wednesday, October 12, 2016 at 6:42:55 PM UTC+2, bartekltg wrote:
>> Nie ma znaczenia.
>> Mnożysz mantysy, które zawsze są w przedziale [0.5,1)
>> cechy dodajesz stałoprzecinkowo.
>> ["Ty w sensie komputer jak mnożysz zmienne float/double",
>> nie trzeba nic ręcznie poprawiać].
>
> Fajny sposób.

Po to był nawias kwadratwy.
To nie jest sposób.
Tak komputer po prostu mnoży liczby zmiennoprzecinkowe.

pzdr
bartekltg

do góry

On Wednesday, October 12, 2016 at 8:15:10 PM UTC+2, bartekltg wrote:
> On 12.10.2016 19:55, M.M. wrote:
> > On Wednesday, October 12, 2016 at 6:42:55 PM UTC+2, bartekltg wrote:
> >> Nie ma znaczenia.
> >> Mnożysz mantysy, które zawsze są w przedziale [0.5,1)
> >> cechy dodajesz stałoprzecinkowo.
> >> ["Ty w sensie komputer jak mnożysz zmienne float/double",
> >> nie trzeba nic ręcznie poprawiać].
> >
> > Fajny sposób.
>
> Po to był nawias kwadratwy.
> To nie jest sposób.
> Tak komputer po prostu mnoży liczby zmiennoprzecinkowe.

No tak, ale w C++ (chyba) nie ma instrukcji która wymnoży w ten sposób
więcej niż dwie liczby? Powiedzmy że mamy 10 dużych liczb i 10 małych.

Gdy zacznę mnożyć od małych, to pewnie najpierw osiągnę zero, więc
potem mnożenie przez duże liczby da wynik zero.

Gdy zacznę mnożyć od dużych, to dojdzie do przepełnienia typu,
pewnie będzie +inf, i mnożenie z małymi zakończy się błędem.

Dobrze byłoby posortować: raz duża, raz mała, ale to trudne jest.

Lepiej wyciągnąć po każdym mnożeniu wykładnik sumować w osobnej
zmiennej, a w bieżącym wyniku go zerować. Coś mniej/więcej tak:

wynik = 1;
E = 0;
tab[N];
for( i=0 ; i<N ; i++ ) {
wynik *= tab[i];
E += wykladnik( wynik );
wynik /= baza ^ wykladnik( wynik );
}
return wynik *= baza ^ E;

Pozdrawiam

do góry

On 12.10.2016 20:59, M.M. wrote:
> On Wednesday, October 12, 2016 at 8:15:10 PM UTC+2, bartekltg wrote:
>> On 12.10.2016 19:55, M.M. wrote:
>>> On Wednesday, October 12, 2016 at 6:42:55 PM UTC+2, bartekltg wrote:
>>>> Nie ma znaczenia.
>>>> Mnożysz mantysy, które zawsze są w przedziale [0.5,1)
>>>> cechy dodajesz stałoprzecinkowo.
>>>> ["Ty w sensie komputer jak mnożysz zmienne float/double",
>>>> nie trzeba nic ręcznie poprawiać].
>>>
>>> Fajny sposób.
>>
>> Po to był nawias kwadratwy.
>> To nie jest sposób.
>> Tak komputer po prostu mnoży liczby zmiennoprzecinkowe.
>
> No tak, ale w C++ (chyba) nie ma instrukcji która wymnoży w ten sposób
> więcej niż dwie liczby? Powiedzmy że mamy 10 dużych liczb i 10 małych.

A ile liczb chcesz mnożyć?

Przeczytaj, jak są zbudowane liczby zmiennoprzecinkowe.

Teraz widzisz?

Operacja pomnożenia ich po kolei (jedna po drugiej)
robi to co opisałem.
Mnoży (też po kolei!) mantysy, a cechy się dodają.


> Gdy zacznę mnożyć od małych, to pewnie najpierw osiągnę zero, więc
> potem mnożenie przez duże liczby da wynik zero.

Spodziewasz się przekroczyć 10^-308?
Wtedy dopiero wchodzimy w zakres liczb podnormalnych
i tracimy.
Jak rzeczywisćie masz taki wtedny szereg, to posortuj
i jeśli dotychczsowy wynik jest <1, bierz liczbę
z "małego" końca, jeśli > 1, bierz z dużego.

Wtedy, jeśli wynik mieści się w zakresie,
też nie powinieneś tego zakresu przekroczyć.


> Dobrze byłoby posortować: raz duża, raz mała, ale to trudne jest.

Wystarczy brać z obu końców ;->

Ale sortowanie jest kosztowne;-)

Jeśli liczby nie są za duże/za małe, wyciagaj co jakiś
czas cechę i trzymaj w jakimś int64_t ;-)

Wyciaga się tym:
http://en.cppreference.com/w/cpp/numeric/math/frexp


> Lepiej wyciągnąć po każdym mnożeniu wykładnik sumować w osobnej
> zmiennej, a w bieżącym wyniku go zerować. Coś mniej/więcej tak:
>
> wynik = 1;
> E = 0;
> tab[N];
> for( i=0 ; i<N ; i++ ) {
> wynik *= tab[i];
> E += wykladnik( wynik );
> wynik /= baza ^ wykladnik( wynik );
> }

1. Używaj frexp.



2. Naprawdę masz aż tak wielkie liczby? I jednocześnie umiarkowany
wynik? Doble trzyma od 10^-308 do 10^308 z taką samą precyzją.

Jeśli liczby mają wykłądniki [-10, 10] spokojnie możesz mnożyć
po 30 w kupie.

[BTW, jeśli to dwumian Newtona, to daje się go rozsądniej obliczać;)]



Wracając do błędów. Co prawda mnożenie jest porządne,
w przeciwieństwie do dodawania nigdy nie tracimy cyfr
znaczacych (uwarunkowanie jest równe 1), to mozę jednak
występuje efekt podobny jak widoczny w porównaniu
normalnego sumowania i sumowania 'parami'.
https://en.wikipedia.org/wiki/Pairwise_summation
Można by sprawdzić, czy dla mnożenia też daje lepszy wynik.

Z linka zerknij na odnośnik [1]. To chyba o tej pracy wspominałem
poprzednio.

Prosty test (kod na końcu)

Wyniki... niejednoznaczne;-)


-2.91802e-13
-3.21668e-12
3.16062e-14
-5.34369e-14

-2.91802e-13
-3.21668e-12
1.26974e-13
-2.03522e-14

Mnożenie parami czasem coś poprawia, czasem nie,
Pogarsza znacznie dla ciągu po kolie.
Za to losowa kolejność poprawia wynik;-)

Elementami ciągu są liczby niewymeirne bliskie jedynce,
iloczyn powinien być równy 1, z powodu błędów zaokrągleń
na poziomie wyliczania wartośći 'dokaladną' sumę licze
zmienną GMP sporej precyzji.

w sumie 10^6 elementów, każdy z błędem rządu 10^-16,
10^-13 ma prawo odparować.


pzdr
bartekltg


#include <iostream>
#include <random>
#include <boost/multiprecision/gmp.hpp>
#include <chrono>
#include <random>
#include <algorithm>

template <typename iterator>
double recmul(const iterator & first, const iterator & second)
{
if (second-first == 1 ) return *first;
else
{
auto mid = first + (second-first)/2;
return recmul(first, mid)*recmul(mid, second);
}


}


int main()
{
const int N = 1024*1024;
vector <double> tab;
tab.reserve(2*N);

for (int i=1; i<N;i++)
{
tab.push_back(sqrt( i/(double)(i+1) ));
tab.push_back(sqrt( (i+1)/(double)i ));
}

boost::multiprecision::mpf_float_1000 MP_wyn=1;
for (auto it = begin(tab); it<end(tab); ++it )
MP_wyn*=*it;



double aku=1.0;
for (auto it = begin(tab); it<end(tab); ++it )
aku*=*it;

cout<< aku-MP_wyn << endl;
cout<< recmul(begin(tab),end(tab))-MP_wyn << endl;

random_device rd;
default_random_engine gen(rd());
shuffle(begin(tab),end(tab),gen);

aku=1.0;
for (auto it = begin(tab); it<end(tab); ++it )
aku*=*it;


cout<< aku-MP_wyn << endl;
cout<< recmul(begin(tab),end(tab))-MP_wyn << endl;

return 0;

}

do góry

On Wednesday, October 12, 2016 at 10:19:04 PM UTC+2, bartekltg wrote:
> On 12.10.2016 20:59, M.M. wrote:
> > On Wednesday, October 12, 2016 at 8:15:10 PM UTC+2, bartekltg wrote:
> >> On 12.10.2016 19:55, M.M. wrote:
> >>> On Wednesday, October 12, 2016 at 6:42:55 PM UTC+2, bartekltg wrote:
> >>>> Nie ma znaczenia.
> >>>> Mnożysz mantysy, które zawsze są w przedziale [0.5,1)
> >>>> cechy dodajesz stałoprzecinkowo.
> >>>> ["Ty w sensie komputer jak mnożysz zmienne float/double",
> >>>> nie trzeba nic ręcznie poprawiać].
> >>>
> >>> Fajny sposób.
> >>
> >> Po to był nawias kwadratwy.
> >> To nie jest sposób.
> >> Tak komputer po prostu mnoży liczby zmiennoprzecinkowe.
> >
> > No tak, ale w C++ (chyba) nie ma instrukcji która wymnoży w ten sposób
> > więcej niż dwie liczby? Powiedzmy że mamy 10 dużych liczb i 10 małych.
>
> A ile liczb chcesz mnożyć?
>
> Przeczytaj, jak są zbudowane liczby zmiennoprzecinkowe.
>
> Teraz widzisz?
>
> Operacja pomnożenia ich po kolei (jedna po drugiej)
> robi to co opisałem.
> Mnoży (też po kolei!) mantysy, a cechy się dodają.
>
>
> > Gdy zacznę mnożyć od małych, to pewnie najpierw osiągnę zero, więc
> > potem mnożenie przez duże liczby da wynik zero.
>
> Spodziewasz się przekroczyć 10^-308?
> Wtedy dopiero wchodzimy w zakres liczb podnormalnych
> i tracimy.
> Jak rzeczywisćie masz taki wtedny szereg, to posortuj
> i jeśli dotychczsowy wynik jest <1, bierz liczbę
> z "małego" końca, jeśli > 1, bierz z dużego.
>
> Wtedy, jeśli wynik mieści się w zakresie,
> też nie powinieneś tego zakresu przekroczyć.
>
>
> > Dobrze byłoby posortować: raz duża, raz mała, ale to trudne jest.
>
> Wystarczy brać z obu końców ;->
>
> Ale sortowanie jest kosztowne;-)
>
> Jeśli liczby nie są za duże/za małe, wyciagaj co jakiś
> czas cechę i trzymaj w jakimś int64_t ;-)
>
> Wyciaga się tym:
> http://en.cppreference.com/w/cpp/numeric/math/frexp
>
>
> > Lepiej wyciągnąć po każdym mnożeniu wykładnik sumować w osobnej
> > zmiennej, a w bieżącym wyniku go zerować. Coś mniej/więcej tak:
> >
> > wynik = 1;
> > E = 0;
> > tab[N];
> > for( i=0 ; i<N ; i++ ) {
> > wynik *= tab[i];
> > E += wykladnik( wynik );
> > wynik /= baza ^ wykladnik( wynik );
> > }
>
> 1. Używaj frexp.
>
>
>
> 2. Naprawdę masz aż tak wielkie liczby? I jednocześnie umiarkowany
> wynik? Doble trzyma od 10^-308 do 10^308 z taką samą precyzją.
>
> Jeśli liczby mają wykłądniki [-10, 10] spokojnie możesz mnożyć
> po 30 w kupie.
>
> [BTW, jeśli to dwumian Newtona, to daje się go rozsądniej obliczać;)]
>
>
>
> Wracając do błędów. Co prawda mnożenie jest porządne,
> w przeciwieństwie do dodawania nigdy nie tracimy cyfr
> znaczacych (uwarunkowanie jest równe 1), to mozę jednak
> występuje efekt podobny jak widoczny w porównaniu
> normalnego sumowania i sumowania 'parami'.
> https://en.wikipedia.org/wiki/Pairwise_summation
> Można by sprawdzić, czy dla mnożenia też daje lepszy wynik.
>
> Z linka zerknij na odnośnik [1]. To chyba o tej pracy wspominałem
> poprzednio.
>
> Prosty test (kod na końcu)
>
> Wyniki... niejednoznaczne;-)
>
>
> -2.91802e-13
> -3.21668e-12
> 3.16062e-14
> -5.34369e-14
>
> -2.91802e-13
> -3.21668e-12
> 1.26974e-13
> -2.03522e-14
>
> Mnożenie parami czasem coś poprawia, czasem nie,
> Pogarsza znacznie dla ciągu po kolie.
> Za to losowa kolejność poprawia wynik;-)
>
> Elementami ciągu są liczby niewymeirne bliskie jedynce,
> iloczyn powinien być równy 1, z powodu błędów zaokrągleń
> na poziomie wyliczania wartośći 'dokaladną' sumę licze
> zmienną GMP sporej precyzji.
>
> w sumie 10^6 elementów, każdy z błędem rządu 10^-16,
> 10^-13 ma prawo odparować.
>
>
> pzdr
> bartekltg
>
>
> #include <iostream>
> #include <random>
> #include <boost/multiprecision/gmp.hpp>
> #include <chrono>
> #include <random>
> #include <algorithm>
>
> template <typename iterator>
> double recmul(const iterator & first, const iterator & second)
> {
> if (second-first == 1 ) return *first;
> else
> {
> auto mid = first + (second-first)/2;
> return recmul(first, mid)*recmul(mid, second);
> }
>
>
> }
>
>
> int main()
> {
> const int N = 1024*1024;
> vector <double> tab;
> tab.reserve(2*N);
>
> for (int i=1; i<N;i++)
> {
> tab.push_back(sqrt( i/(double)(i+1) ));
> tab.push_back(sqrt( (i+1)/(double)i ));
> }
>
> boost::multiprecision::mpf_float_1000 MP_wyn=1;
> for (auto it = begin(tab); it<end(tab); ++it )
> MP_wyn*=*it;
>
>
>
> double aku=1.0;
> for (auto it = begin(tab); it<end(tab); ++it )
> aku*=*it;
>
> cout<< aku-MP_wyn << endl;
> cout<< recmul(begin(tab),end(tab))-MP_wyn << endl;
>
> random_device rd;
> default_random_engine gen(rd());
> shuffle(begin(tab),end(tab),gen);
>
> aku=1.0;
> for (auto it = begin(tab); it<end(tab); ++it )
> aku*=*it;
>
>
> cout<< aku-MP_wyn << endl;
> cout<< recmul(begin(tab),end(tab))-MP_wyn << endl;
>
> return 0;
>
> }

Dzięki za odpowiedź, najbardziej przyda się pomysł
który mi podsunąłeś, żeby wyciągać eksponenty :)
Resztę w miarę ogarniam.

Mam takie liczby do wymnożenia, że przy 80 losowo wybranych
często pada. Liczby, choć nie służą do wyliczenia dwumianiu
newtona, to powstają z dużych dodatnich i ujemnych potęg, więc
mogą po wymnożeniu nawet dać coś w okolicach jedynki. Napisałem
najpierw podobnie do sortowania. Czynniki mam w kilku wektorach,
jeśli w jednym mam czynniki duże, to w drugim mam małe, albo
na odwrót. S skrócie: próbuję wymnożyć przez jedną, potem przez
drugą i wybieram tę, która daje wynik bliższy jedynce. Po tym
zabiegu mogę wymnożyć nawet 400 liczb bez problemów z zakresem, a
nie muszę sortować. Jednak pomysł z przeniesieniem wykładników
do osobnej zmiennej zapewne będzie najlepszy i szybki.

Pozdrawiam

do góry

On Thursday, October 13, 2016 at 3:41:25 PM UTC+2, M.M. wrote:

>
> Dzięki za odpowiedź, najbardziej przyda się pomysł
> który mi podsunąłeś, żeby wyciągać eksponenty :)
> Resztę w miarę ogarniam.



> Mam takie liczby do wymnożenia, że przy 80 losowo wybranych
> często pada. Liczby, choć nie służą do wyliczenia dwumianiu
> newtona, to powstają z dużych dodatnich i ujemnych potęg, więc
> mogą po wymnożeniu nawet dać coś w okolicach jedynki. Napisałem
> najpierw podobnie do sortowania. Czynniki mam w kilku wektorach,
> jeśli w jednym mam czynniki duże, to w drugim mam małe, albo
> na odwrót. S skrócie: próbuję wymnożyć przez jedną, potem przez
> drugą i wybieram tę, która daje wynik bliższy jedynce. Po tym
> zabiegu mogę wymnożyć nawet 400 liczb bez problemów z zakresem, a
> nie muszę sortować. Jednak pomysł z przeniesieniem wykładników
> do osobnej zmiennej zapewne będzie najlepszy i szybki.

Zrobiłem mały test
http://pastebin.com/3mWyh93z

Mnożenie bezpośrednie (najczęśceij daje 0 lub inf ;-))
Dodanie logartrymów - po to, aby sprawdzić wynik. Niestabilne!
Ale daje wynik i wiadomo, że sie nie pomyliśmy.

Prawdziwe algorytmy:

Jeden to modyfikacja sortowania. Nie potrzeba nam przecież
sortować, wystarczy zrobić podział na > 1 i < 1.
Potem mnożymy, jeśli tymczasowy wynik jest >1, korzystamy
ze zbioru liczb mniejszych, jeśli mniejszy, z większych.
Potem trzeba domnożyć ogon.

Rozwinięciem tego jest zrobienie partition i mnożenia na raz.
Mam dwa iteratory, jeden omija wszystkie <1, drugi omija
wszystko co jest >1.
Dalej tak samo, używam tego z przecinwej strony niz wynik.

Różnica jest minimalna.
frexp jest ponad 2 razy szybszy!

W sumie mogliśmy się tego spodziewać. W obu wersjach 'sortujących'
mamy dwa przebiegi po tablicy, dla frexp tylko jeden,
a to dostęp do pamięci, a nie operacje będą wąskim gardłem.

Chyba, że masz tylko tyle liczb, ile się mieści w cache.
Wtedy... kto wie, trzeba by mój test nieco przerobić.


pzdr
bartekltg

do góry

W dniu środa, 12 października 2016 17:41:35 UTC+2 użytkownik M.M. napisał:
> Mamy dane wejściowe:
> tab[N] - tablica liczb zmiennoprzecinkowych
>
W operacjach można obawiać się albo zaniku wartości, albo przepełnienia. Posortuj
więc rosnąco co wartosci abs, a następnie mnóż największe z nich przez najmniejsze
(co do abs, rzecz jasna)
>
> Pozdrawiam

do góry