Numerička matematika

Današnja računala obično koriste 32-bitni (float) i 64-bitni (double) IEEE floating point prikaz kako bi prikazali realne brojeve. Za više informacija o samom prikazu najbolje je pogledati sljedeće odlične Topcoder tutoriale:

http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=integersReals
http://community.topcoder.com/tc?module=Static&d1=tutorials&d2=integersReals2

Cilj ovog članka je da shvatite probleme koji se mogu javiti oko preciznost kada radite s realnim brojevima. Uvesti ćemo pojam apsolutne i relativne pogreške koji su potpuno centralni za ovu temu. U nekoj realnoj situaciji jedina stvar oko koje ćete se brinuti pri radu s realnim brojevima će biti da vam apsolutna ili relativna pogreška (što god je manje) bude manja od npr. .

Relativna i apsolutna pogreška

Pretpostavimo da je stvarna vrijednost neke veličine , a mi je u programu reprezentiramo s varijablom čija je stvarna vrijednost . Tada definiramo:

Nekakva normalna vrijednost za min(relativna pogreška, apsolutna pogreška) je reda veličine .

Primjeri

Nekoliko najkorisnijih primjera:

  1. Zbrajanje: .
  2. Oduzimanje: .
    1. Dakle ako , relativna pogreška oduzimanja može eksplodirati.
  3. Množenje: .
  4. Dijeljenje: .
Izvod

Pretpostavimo u ovom trenutku kako su računala savršeno precizna i da su sve pogreške unesene mjerenjima i aproksimacijama. Zanima nas sljedeći problem: za zadane brojeve izračunati za neku "normalnu" (barem dva puta derivabilnu) funkciju .

Ovdje smo zanemarili članove višeg reda. Ovo je opravdano jer su oni normalno reda veličine MathJax error, mnogo manji nego što je reprezentabilno. Sada imamo:

Zadatak: CollectingBonuses

Izvor: TopCoder SRM 400; Div I; Level 3
Zadatak: http://community.topcoder.com/stat?c=problem_statement&pm=8764&rd=12172&rm=&cr=10574855
Editorial: http://community.topcoder.com/tc?module=Static&d1=match_editorials&d2=srm400

Nakon što malo razmislite o zadatku (i znate ponešto o Markovljevim lancima), nije teško zaključiti kako je potrebno izračunati sljedeću sumu:

Problem nam tvore ograničenja: . Naivno sumiranje for-petljom ne dolazi u obzir. Problemu je potrebno drugačije pristupiti: definirajmo takozvani harmonijski red:

Sada je potrebno izračunati s relativnom ili apsolutnom pogreškom od (dovoljno je ispuniti bilo koji od ta dva uvjeta). Kratko prisječanje na gornje formule nam govori da je to zadovoljeno kada izračunamo razliku harmonijskih suma s relativnom pogreškom od ili s apsolutnom pogreškom od (dovoljno je ispuniti bilo koji od ta dva uvjeta). No, broj s apsolutnom pogreškom od nije nemoguće precizno pohraniti (zato jer harmonijske sume nisu brojevi bliski 0).

Računanje harmonijskog reda je poznat problem u matematici i njemu se pristupa tako da se usporedi s integralom u što nećemo ulaziti. WolframAlpha nam je ovdje od velike pomoći i kaže da:

gdje je . Iako ova formula ima veliku relativnu i apsolutnu pogrešku za male , primjetite da sada uvijek znamo izračunati s malom relativnom pogreškom. Kako je za sve praktične primjene manji od 100, dovoljno je izračunati s aps. pogreškom :

ako je n < 10^3, iskoristi naivni algoritam
inače vrati $\ln n + \gamma + 1/(2n) - 1/(12n^2) + 1/(120n^4)$

No, to NIJE dovoljno za riješiti cijeli zadatak. Budući da oduzimamo dvije harmonijske sume (koje su obično bliske za velike n, k), relativna pogreška se može jako povećati, stoga to nije konačan način na koji se zadatak rješava. U tome slučaju kada su oba n i k jako veliki, a njihove harmonijske sume bliske, koristiti ćemo formulu:

Još jednom, ova formula za razliku je neprecizna kada je dovoljno mali broj (aproksimacija harmonijske sume je neprecizna za male H). Ja (Goran) sam osobno ručno izračunao prvih milijun članova sume, pa radio aproksimacije za preostale članove (naknadno sam otkrio da gornji algoritam prolazi čak i kada se ručno posumiraju samo prvih 10 članova).

Potrebno je dati još samo jednu napomenu: funkcija , koja je u C/C++-u implementirana u math.h preko funkcije double log(double) neprecizna ukoliko nam je argument funkcije vrlo blizak 1 (npr. ). U tome slučaju je potrebno koristiti funkciju log1p(x) koja nam vraća , ali s mnogo većom preciznosti.

U konačnici, algoritam izgleda ovako (nisam ga detaljno provjerio, ali trebao bi biti ok):

ako k < 10^6 vrati n*naivno_sumiranje(n-k, n);
prefix_suma := naivno_sumiranje(n-k+1, n-k+10^6)  // ručno posumiraj prvih 10^6 članova
k := k - 10^6
ako n/(n-k) < 1.1 tada L := log1p(k/(n-k)) inače L := log(n/(n-k))
vrati n * (prefix_suma + L - k/2/n/(n-k) + k*n/12/n^2/(n-k)^2)
© 2012. Goran Žužić Creative Commons License

Ovaj članak objavljen je pod
Creative Commons Attribution-ShareAlike 3.0 Croatia License