Re: Moninkertaisen tarkkuuden laskentaa

[vastaus aiempaan viestiin]

Saamatta toistaiseksi juuri palautetta muilta kuin itseltäni :)
olen edelleen jatkanut puuhailua "suurten lukujen" parissa.
ARIT-operaatio kykenee nyt laskemaan myös standardifunktioiden
exp,log,sin,cos,tan ja arctan arvoja mielivaltaisella tarkkuudella:
.....................................................................
Esim.  ACCURACY=50
ARIT [CUR+1]=arctan(10)
0.14711276743037345918528755717617308518553063771832e1
ARIT [CUR+1]=tan([CUR-1])
0.99999999999999999999999999999999999999999999999961e1
.....................................................................
Myös editoriaalinen laskenta (melkein kaikilta osin) toimii monin-
kertaisella tarkkuudella käyttämällä ACCURACY-arvoa, joka on yli 16.
.....................................................................
Esim. ACCURACY=100
1/3=0.3333333333333333333333333333333333333333333333333333333333333333\
333333333333333333333333333333333333
potenssi=2^n     n=100
potenssi=1267650600228229401496703205376
potenssi^(1/n)=2
.......................................................................
Katsotaan miten (1+1/n)^n lähestyy suurella n-arvolla lukua e=exp(1):
n=10^500  ACCURACY=1000,10  (10 tarkoittaa tulostustarkkuutta)
exp(1)-(1+1/n)^n=0.135914091e-499
.......................................................................
Geometrisen sarjan summa 1+1/2+1/4+1/8+1/16+... = 2
ACCURACY=1000,50
summa=for(i=1)to(2000)term(T=1)sum(0.5*T)
  2-summa=0.1741961963243443335115239098955774459171820748541e-601
Tarkistus laskemalla jäännöstermi:
2^(-1999)=0.1741961963243443335115239098955774459171820748541e-601
.......................................................................

Laskenta suurilla tarkkuuksilla näyttää olevan oma "taiteenlajinsa".
Tehokkaaseen toimintaan tarvitaan joskus aika erikoisia temppuja.
Hyvä esimerkki on logaritmien laskenta, joka (vaadittaessa todella
suurta tarkkuutta) näyttää parhaiten onnistuvan seuraavasti:
.......................................................................
Luvun 10 luonnollinen logaritmi tarkkuudella ACCURACY=50 on
log(10)=0.23025850929940456840179914546843642076011014886287e1
SURVO MM laskee sen itse asiassa seuraavalla kaavalla:
L(x,m):=#PI/(2*m*agm(1,4/x^m))
josta saadaan tuloksena
L(10,26)=0.23025850929940456840179914546843642076011014886287e1
.......................................................................
Tässä #PI on tietenkin tuo tunnettu vakio 3.14159265... laskettuna
ja talletettuna pysyvästi tiedostoon
<Survo>\U\SYS\#PI.NBR ainakin 100000 merkitsevällä numerolla.
Se on tehty ARIT-komennolla
ARIT PI,#PI / ACCURACY=100000 ja mukana Survo-asennuksessa.
Parametrin m on oltava suurempi kuin ACCURACY/2.
agm(x,y) on lukujen x ja y aritmeettis-geometrinen keskiarvo ja se
lasketaan seuraavalla iterointimenettelyllä:

Olkoon x(0)=x, y(0)=y.
Iterointiaske vaiheesta n vaiheeseen n+1:
x(n+1)=[x(n)+y(n)]/2       aritmeettinen keskiarvo
y(n+1)=sqrt[x(n)*y(n)]     geometrinen keskiarvo
Iterointia jatketaan, kunnes x(n)=y(n).

AGM-tekniikalla saadaan myös pi:n likiarvoja yllättävän nopeasti.
Tästä löytyy esimerkki PITKARIT-toimituskentästä, joka tulee olemaan
mukana SURVO MM-asennuksessa.
Yleensä standardifunktiot (kuten exp,sin,...) kuitenkin lasketaan
sarjakehitelmistä, joissa suppenemista terästetään sopivilla
argumentin muunnoksilla.

Kuten jo kuukauden takaisessa viestissäni totesin, tavallinen
kaksoistarkkuus (15-16 merkitsevää numeroa) riittää mainiosti
esim. tyypillisiin tilastollisiin laskelmiin ja analyyseihin.
ARIT-operaation ja moninkertaisen tarkkuuden editoriaalisen
laskennan avulla on (kuten jo edellä olevista esimerkeistä ilmenee)
mahdollista havainnollistaa erilaisia numeerisia prosesseja
kuten sarjojen suppenemista ja lausekkeiden lähestymistä kohti
raja-arvojaan.

On myös aitoja ongelmia, joissa tarvitaan parempaa laskutarkkuutta ja
mahdollisuutta käsitellä yli 10^308 ja alle 10^-308 suuruusluokkaa
olevia lukuja.
Näytteenä tällaisesta käsittelen seuraavaa leipurin ongelmaa:

                     *********

Taikinasta leivotaan 1000 rusinapullaa. Siihen on sekoitettuna 2000
rusinaa. Mikä on todennäköisyys, että jokaiseen pullaan tulee ainakin
yksi rusina?

Ennenkuin luet eteenpäin, arvioi mielessäsi, mitä suuruusluokkaa ko.
todennäköisyys voisi olla? Onko se edes sadasosa?

Ratkaisu:

Rusinapullia N kpl. (1000)
Rusinoita M kpl. (2000)

Määritellään tapahtumat
Ai={Pullassa i ei rusinoita} i=1,2,...,N

1-P(A1+A2+...+AN) on tn. että jokaisessa pullassa ainakin yksi rusina.

Koska tapahtumat A1,A2,...,AN eivät ole toisensa poissulkevia, on
käytettävä todennäköisyyslaskennan yleistä yhteenlaskulausetta
P(A1+A2+...+AN) = P(A1)+P(A2)+...+P(AN) -sum P(AiAj) +sum P(AiAjAk) -...

Tässä tapauksessa saadaan
P(Ai)=[(N-1)/N]^M,  i=1,2,...,N
P(AiAj)=[(N-2)/N]^M kaikille i <> j
P(AiAjAk)=[(N-3)/N]^M kaikille i <> j <> k

ja siis
P(A1+A2+...+AN) =
    N*[(N-1)/N]^M -C(N,2)*[(N-2)/N]^M +C(N,3)*[(N-3)/N]^M -...
.......................................................................
Jos näitä todennäköisyyksiä yrittää laskea tavallisella editoriaalisella
laskennalla, haluttu todennäköisyys 1-P(N,M) on määrättävissä
suurilla rusinamäärillä seuraavasti:

P(N,M):=for(i=1)to(N-1)sum(S(N,M,i)*(-1)^(i+1))
S(N,M,i):=exp(lfact(N)-lfact(i)-lfact(N-i)+M*(log(N-i)-log(N)))
Tässä binomikerroin C(N,i) etc. joudutaan laskemaan logaritmien kautta,
jotta vältyttäisiin ylisuurilta luvuilta.
Esim.
1-P(1000,9874)=0.95003468545475
mikä on oikein 9 desimaalin tarkkuudella.
(Rusinoita tarvitaan siis lähes kymmenen tuhatta, jotta jokaiseen
pullaan tulisi edes yksi rusina todennäköisyydellä 0.95)

Nyt kuitenkin kysytyssä tapauksessa M=2000 saataisiin
1-P(1000,2000)=1.0980563252807e+037
eli aivan mieletön tulos ("todennäköisyys" yli 10^37 !!!).
Tämä johtuu siitä, että summa koostuu vaihtuvanmerkkisistä termeistä,
jotka ovat pienillä M-arvoilla itseisarvoltaan erittäin suuria.
Pienellä tarkkuudella ne kumoavat toisiaan niin, kaikki merkitsevät
numerot menetetään.
.......................................................................
Sama tehdään nyt tarkkuudella ACCURACY=200,15 .

P(N,M):=for(i=1)to(N-1)sum(S(N,M,i)*(-1)^(i+1))
S(N,M,i):=C(N,i)*((N-i)/N)^M

1-P(1000,9874)=0.95003468545506
1-P(1000,2000)=0.18643919830164e-77
eli kysytyssä 2000 rusinan tilanteessa ko. todennäköisyys on
huikean pieni (siis suuruusluokkaa 10^-78).
.......................................................................
Vielä hankalammaksi laskenta tulee pienimmällä mahdollisella arvolla
M=1000.

P(N,M):=for(i=1)to(N-1)sum(S(N,M,i)*(-1)^(i+1))
S(N,M,i):=C(N,i)*((N-i)/N)^M

ACCURACY=600,15 / Tässä joudutaan käyttämään vielä suurempaa tarkkuutta.
TIME COUNT START
1-P(1000,1000)=0.40238726007709e-432
TIME COUNT END   345.507

ja aikaakin kului lähes 6 minuuttia 1.6 GHz:n suorittimella.

Todennäköisyys sille, että 1000 pullasta jokaiseen tulee rusina, kun
rusinoita on käytettävissä vai yksi kutakin pullaa kohti, on siis
suuruusluokkaa 10^-432.

Tässä erikoistapauksessa tulos on helposti tarkastettavissa, sillä
nyt ko. todennäköisyys voidaan esittää yksinkertaisesti muodossa

(N-1)/N*(N-2)*N*(N-3)/N*...*1/N = N!/N^N eli

fact(1000)/1000^1000=0.40238726007709e-432

.......................................................................
Kysytylle todennäköisyydelle on kehitetty myös likiarvokaava
P2(N,M):=1-exp(-N*exp(-M/N))
mutta se toimii vain suurilla M-arvoilla tyydyttävästi:
ACCURACY=600,15
1-P2(1000,9874)=0.94980713501579       (oikea suuruusluokka)
1-P2(1000,1000)=0.17060379734841e-159  (virheellinen suuruusluokka)
.......................................................................

                     *********

Sekä suuren tarkkuuden editoriaalinen laskenta että ARIT-operaatiot
tulevat olemaan mukana SURVO MM:n versioista 1.24 lähtien.
Vaikka olen simputtanut näitä toimintoja monilla tavoilla ja
saanut korjatuksi lukuisia virheitä, on toistaiseksi oltava
varovainen niitä käytettäessä.
Toivon, että käyttäjät kertoisivat välittömästi havaitsemistaan
ongelmista.
Olen verrannut näiden toimintojen nopeutta eräiden matematiikka-
ohjelmien (Maple V ja Derive) vastaaviin. Varsinkin todella suurilla
tarkkuuksilla SURVO MM on joissain toiminnoissa (esim. sqrt ja log)
näitä nopeampi, joissain (esim. exp ja arctan) taas hitaampi.

Käyttämäni algoritmit perustuvat mm. seuraaviin lähteisiin:

David H. Bailey: A Portable High Performance Multiprecision Package.
RNR Technical Report RNR-90-022 (1993)

Richard P. Brent: Fast Multiple-Precision Evaluation of Elementary
Functions. Journal of the ACM, Vol.23, No.2, April 1976, pp.242-251.

Donald E. Knuth: The Art of Computer Programming, Vol.2.
Addison Wesley, 1981.

Luettelo keskustelunaiheista

Survo-keskustelupalstan (2001-2013) viestit arkistoitiin aika ajoin sukrolla, joka automaattisesti rakensi viesteistä (yli 1600 kpl) HTML-muotoisen sivukokonaisuuden. Vuoden 2013 alusta Survo-keskustelua on jatkettu entistäkin aktiivisemmin osoitteessa forum.survo.fi. Tervetuloa mukaan!