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;

}