Autor Wątek: Liczby Bernoulliego - algorytm  (Przeczytany 751 razy)

Offline beem

  • Użytkownik
    • Wordpress Blog

# Styczeń 09, 2017, 13:12:37
Mam problem z uzyskaniem kolejnych liczb Bernoulliego. Używam algorytmu Akiyama–Tanigawy, który u mnie wygląda tak:
double bernoulli_number(unsigned int n)
{
double* A = new double[n+1];

for (unsigned int m = 0; m <= n; m++)
{
A[m] = 1.0 / (m + 1.0);

for (unsigned int j = m; j > 0; j--)
{
A[j - 1] = j*(A[j - 1] - A[j]);
}
}

double res = A[0];

delete[] A;

return res;
}

Algorytm działa, ale każda kolejna liczba zaczyna być opatrzona coraz to większym błędem. Np dla 30 kolejnych liczb otrzymuje się:
1 0.5
2 0.1666666666666666
3 -2.220446049250313e-16
4 -0.03333333333333399
5 1.77635683940025e-15
6 0.02380952380956758
7 3.794742298168785e-13
8 -0.03333333333154209
9 -1.044986319698182e-11
10 0.07575757526731397
11 -9.93349424938117e-09
12 -0.2531137033204133
13 -1.556503844035717e-06
14 1.166667679153587
15 0.0006419491273961242
16 -7.068502074522296
17 0.6405601281966771
18 69.99698015616998
19 323.0140711709549
20 5986.641669860796
21 120040.0651014931
22 1589860.901390899
23 -19576696.20832111
24 -3163545259.292155
25 -195508362222.1589
26 -9583584454807.713
27 -413822699863629.6
28 -1.630954724494887e+16
29 -5.98305130690999e+17
W miejscach nieparzystych 3, 5, 7 itd powinno być 0.0. Np B18 = 54.97 a u mnie 69.99.
Ktoś wie co może być przyczyną tak szybko rosnącego błędu? Wydaje mi się, że to jest błąd zaokrąglenia, ale może są to jakieś ustawienia kompilatora.
« Ostatnia zmiana: Styczeń 09, 2017, 20:55:34 wysłana przez beem »

Offline Mr. Spam

  • Miłośnik przetworów mięsnych

Offline mwojt

  • Użytkownik

# Styczeń 09, 2017, 17:30:27
Widzę że pojawiają się bardzo małe liczby: e-16 więc pewnie gdzieś kończy się precyzja double.

Offline Adam27

  • Użytkownik

# Styczeń 10, 2017, 13:02:29
Precyzja double kończy się już na samym początku. Przy obliczaniu A[3] program wykonuje działanie 1/6 - 2*(1/3-1/4), które nie oblicza się do 0 nawet zakładając brak błędów w poprzednich działaniach: http://ideone.com/KP0dn0. Nic z tym nie zrobisz, bo 1/3 nie ma dokładnej reprezentacji w systemie binarnym.

Lepiej byłoby Ci operować na liczbach wymiernych, np. z użyciem https://gmplib.org/, albo samemu pisząc proste sprowadzanie do wspólnego mianownika (co jednak może z drugiej strony szybko wyczerpać precyzję intów).

Offline beem

  • Użytkownik
    • Wordpress Blog

# Styczeń 10, 2017, 21:17:13
Wygląda to rzeczywiście na małą precyzję double. Winne jest tu dzielenie w pierwszej pętli for. Na szczęście już sobie z tym problemem poradziłem, użyłem innego algorytmu. GMP nie działa z chyba Visualem a szkoda.