Moninkertaisen tarkkuuden laskentaa

[viesti Survo-keskustelupalstalla (2001-2013)]

Kirjoittaja: Seppo Mustonen
Sähköposti:    -
Päiväys: 21.7.2002 18:05

Survossa numeerinen laskenta tapahtuu yleensä liukuluvuilla
noin 15-16 merkitsevällä numerolla eli ns. kaksoistarkkuudella.
Pitemmissä laskelmissa syntyvät tulokset eivät välttämättä ole
näin tarkkoja johtuen matkan varrella tapahtuneista pyöristyksistä.
Tämä laskentatarkkuus on kuitenkin täysin riittävä "kaikkiin" järkeviin
tehtäviin (esim. tilastollisissa sovelluksissa).

Oma viehätyksensä, jota en halua kiistää, on kuitenkin siinä, että
tarvittaessa voitaisiin tehdä laskelmia suuremmallakin tarkkuudella.
Tässä suhteessa ns. matematiikkaohjelmat (esim. Maple ja Mathematica)
ovat varsin eteviä.

Ajatellen erityisesti opetustilanteita päätin eräänlaisena
kesälomaharrastuksena laatia pienen lisäyksen Survoon, joka
sallii peruslaskutoimitukset "mielivaltaisella" tarkkuudella.
On täysin mahdollista ujuttaa tämä laajennus suoraan editoriaaliseen
laskentaan (asetetaan esim. ACCURACY=1000), mutta toistaiseksi olen
tyytynyt toteutukseen erillisen ARIT-komennon muodossa.

Esim. toimituskentän riveille A ja B kirjoitettujen lukujen tulo
saadaan riville C seuraavalla ARIT-komennolla:

*
*ARIT [C]=[A]*[B]
*
A9999999999999999999999
*
B5555555555555555555555
*
C55555555555555555555544444444444444444444445
*

Lukuihin viitataan siis hakasuluissa olevin rivinumeroin tai -tunnuksin
(esim. [35],[A],[B],[CUR-2],[CUR+1],[END+2]).
....................................................................
(Lyhyitä) lukuja voi kirjoittaa myös suoraan komennon sisään.
Esim. 2^2048 (2 potenssiin 2048) lasketaan näin:

*ARIT [CUR+1]=2^2048
*0.32317006071311007300714876688669951960444102669715e617

Oletustarkkuus on 50 merkitsevää numeroa ja tässä tulos tarkoittaa
liukulukua 0.32317006...*10^617 .
.....................................................................
Jos tulos haluttaisiin saada täydellä tarkkuudella, on käytettävä
ACCURACY-täsmennystä (edellisen tuloksen perusteella) arvolla, joka
on ainakin 617 eli siis

*ARIT [X]=2^2048  / ACCURACY=617  WIDTH=50 (tulostusleveys)
X32317006071311007300714876688669951960444102669715\   Pitkät
*48403213034542752465513886789089319720141152291346\   luvut
*36887179609218980194941195591504909210950881523864\   lohkoutuvat
*48283120630877367300996091750197750389652106796057\   näin
*63838406756827679221864261975616183809433847617047\   sopivan
*05816458520363050428875758915410658086075523991239\   mittaisiin
*30385521914333389668342420684974786564569494856176\   pätkiin.
*03532632205807780565933102619270846031415025859286\
*41771167259436037184618573575983511523016459044036\
*97613233287231227125684710820209725157101726931323\
*46967854258065669793504599726835299863821552516638\
*94373355436021354332296046453184786049521481935558\
*53611059596230656
*
*Koska kyseessä on "kakkosen potenssi kakkosen potenssiin",
*tulos [X] voidaan "tarkistaa" ottamalla siitä toistuvasti
*neliöjuuria:
*
*ARIT [CUR+1]=sqrt([X])
*17976931348623159077293051907890247336179769789423\
*06572734300811577326758055009631327084773224075360\
*21120113879871393357658789768814416622492847430639\
*47412437776789342486548527630221960124609411945308\
*29520850057688381506823424628814739131105408272371\
*63350510684586298239947245938479716304835356329624\
*224137216
*ARIT [CUR+1]=sqrt([CUR-7])
*13407807929942597099574024998205846127479365820592\
*39337772356144372176403007354697680187429816690342\
*76900318581864860508537538828119465699464336490060\
*84096
*ARIT [CUR+1]=sqrt([CUR-4])
*11579208923731619542357098500868790785326998466564\
*0564039457584007913129639936
*ARIT [CUR+1]=sqrt([CUR-2])
*340282366920938463463374607431768211456
*ARIT [CUR+1]=sqrt([CUR-1])
*18446744073709551616
*ARIT [CUR+1]=sqrt([CUR-1])
*4294967296
*ARIT [CUR+1]=sqrt([CUR-1])
*65536
*ARIT [CUR+1]=sqrt([CUR-1])
*256
*ARIT [CUR+1]=sqrt([CUR-1])
*16
*ARIT [CUR+1]=sqrt([CUR-1])
*4
*ARIT [CUR+1]=sqrt([CUR-1])
*2

eli 11-kertaisella neliöjuurella palataan tarkasti lukuun 2,
joten eksponenttina oli siis 2^11=2048.
.....................................................................
ARIT-operaatiolla voi havainnollista mukavasti erilaisia laskenta-
algoritmeja. Esim. luvun v neliöjuuren käänteisarvo u lasketaan
Newtonin menetelmällä (välttäen raskaita jakolaskuja) toistamalla
(iteroimalla) kaavaa
u(n+1)=0.5*u(n)*(3-v*u(n)*u(n)),  lim u(n)=u=1/sqrt(v).

Seuraava laskentakaavio toteuttaa tehtävän jatkuvalla aktivoinnilla:

*ACCURACY=400
v 2 / juurrettava
*
*
u 1 / Tähän ilmaantuu 1/sqrt([v])
*
*
*
*
*
*
*Aktivoi seuraava rivi napein F2 ESC:  (Keskeytä painamalla '.')
xARIT [a]=[u]*[u]
*ARIT [a]=[v]*[a]
*ARIT [a]=3-[a]
*ARIT [a]=[u]*[a]
*ARIT [u]=0.5*[a]
*GOTO v-6,x        / Hyppy takaisin riville x
a
*
*
*
*

Kaavion jatkuva aktivointi tuottaa tuloksen:

*ACCURACY=400
v 2 / juurrettava
*
*
u0.70710678118654752440084436210484903928483593768847403658833986899536\
*6239231053519425193767163820786367506923115456148512462418027925368606\
*3220607485499679157066113329637527963778999752505763910302857350547799\
*8580298513726729843100736425870932044459930477616461524215435716072541\
*9881301813997625703994843626698273165904414820310307629176197527372875\
*143879980864917787610168765928505677187301704249423
*
*Aktivoi seuraava rivi napein F2 ESC:  (Keskeytä painamalla '.')
xARIT [a]=[u]*[u]
*ARIT [a]=[v]*[a]
*ARIT [a]=3-[a]
*ARIT [a]=[u]*[a]
*ARIT [u]=0.5*[a]
*GOTO v-6,x        / Hyppy takaisin riville x
a0.49999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999991
*
*Tuloksen tarkistus kertolaskulla:
*ARIT [CUR+1]=[u]*[u]
*0.49999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999999999999999999999999\
*9999999999999999999999999999999999999999999999999991
*

u^2 on siis "riittävällä" tarkkuudella 1/v=1/2.

ARIT laskee neliöjuuret juuri tällä algoritmilla käyttäen lähtöarvona
ko. luvun kaksoistarkkudella laskettua neliöjuurta
(eli tässä tapauksessa 1/sqrt(2)=0.70710678118655).

.....................................................................
Luvun pi likiarvojen määrääminen on ollut kiinnostuksen kohteena kautta
vuosisatojen. Tehtävään löytyy nykyisin hämmästyttävän tehokkaita
algoritmeja.
Näistä on ARIT-ohjelmaan liitetty Borweinin esittämä "Algoritmi 1", kts.
www.cecm.sfu.ca/organics/papers/borwein/paper/html/node4.html 
ja esim. 1000 ensimmäistä desimaalia saadaan alta sekunnissa
(yli 1 GHz:n PC:llä) seuraavasti:

*ARIT PI,CUR+1 / ACCURACY=1001
*0.31415926535897932384626433832795028841971693993751058209749445923078\
*1640628620899862803482534211706798214808651328230664709384460955058223\
*1725359408128481117450284102701938521105559644622948954930381964428810\
*9756659334461284756482337867831652712019091456485669234603486104543266\
*4821339360726024914127372458700660631558817488152092096282925409171536\
*4367892590360011330530548820466521384146951941511609433057270365759591\
*9530921861173819326117931051185480744623799627495673518857527248912279\
*3818301194912983367336244065664308602139494639522473719070217986094370\
*2770539217176293176752384674818467669405132000568127145263560827785771\
*3427577896091736371787214684409012249534301465495853710507922796892589\
*2354201995611212902196086403441815981362977477130996051870721134999999\
*8372978049951059731732816096318595024459455346908302642522308253344685\
*0352619311881710100031378387528865875332083814206171776691473035982534\
*9042875546873115956286388235378759375195778185778053217122680661300192\
*78766111959092164201989e1

                            ***

ARIT-ohjelman toteutus ei ole käytännössä tehokkain mahdollinen,
mutta mm. opetustehtäviin täysin riittävä.
Laskutoimitukset tehdään puhtaina koneen kokonaislukuoperaatioina
käyttäen lukua 1000 kantalukuna.
Lukuja hallitaan liukulukumuodossa siten, että suurin mahdollinen
merkitsevien desimaalinumeroiden määrä on noin 10 miljoonaa.
Suurin mahdollinen luku on noin 10^(2^31-3)=10^2147483645 ja
pienin mahdollinen vastaavasti 10^(-2^31-3)=10^-2147483645.

Avainasemassa nopeuden kannalta on erityisesti se, miten kertolaskut
tehdään. Kun kerrottaessa kaksi n-numeroista lukua keskenään, tarvitaan
tavallisella "allekkainkertomisella" n*n kpl. yksittäisiä kertolaskuja
+ tietty määrä yhteenlaskuja, monille lienee yllätys, että suurilla
n-arvoilla tuo työmäärä voidaan kutistaa lähes kertalukuun n*log(n)
kertaluvun n*n asemasta ja tämä on todella huikea parannus.
Se tapahtuu esim. nopean Fourier-muunnoksen (FFT) avulla.
Verkosta löytyy tietoja mm. hakusanoilla "fast multiplication".

ARIT pystyy tällä hetkellä vain peruslaskutoimituksiin +,-,*,/
sekä kokonaislukupotenssien (^), neliöjuuren (sqrt), sen käänteis-
arvon (invsqrt) sekä pi:n laskentaan. Mm. erilaisia muita alkeis-
funktioita voi tietenkin laskea esim. sarjakehitelmistä, mutta se
on hidasta, ellei niitä ohjelmoida suoraan ARIT:in sisään.

Jään odottamaan käyttäjien reaktioita, ennenkuin päätetään, kannattaako
toimintoja vielä laajentaa ja onko syytä ympätä tämä myös
editoriaaliseen laskentaan, jolloin mutkikkaitakin lausekkeita
saatetaan laskea suoraan hajottamatta niitä perustoimituksiin.

Survon editoriaalinen käyttöliittymä on jälleen ollut omiaan
tätäkin toimintaa kehiteltäessä. Uskon, että liittymän joustavuus
tuntuu myös ARIT-sovelluksissa.

ARIT-moduli tulee olemaan mukana SURVO MM:n versiosta 1.23 alkaen.

Seppo Mustonen

Vastaukset:

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!

Etusivu  |  Keskustelu
Copyright © Survo Systems 2001-2013. All rights reserved.
Updated 2013-06-15.