i i “Klavzar-racunanje-vrednosti” — 2010/5/27 — 10:20 — page 1 — #1 i i i i i i List za mlade matematike, fizike, astronome in računalnikarje ISSN 0351-6652 Letnik 14 (1986/1987) Številka 2 Strani 116–124 Sandi Klavžar in Matija Lokar: RAČUNANJE VREDNOSTI NEKATERIH MATEMA- TIČNIH FUNKCIJ Ključne besede: matematika, računalništvo, algoritmi, elementarne funkcije, kalkulator HP. Elektronska verzija: http://www.presek.si/14/826-Lokar.pdf c© 1986 Društvo matematikov, fizikov in astronomov Slovenije c© 2010 DMFA – založništvo Vse pravice pridržane. Razmnoževanje ali reproduciranje celote ali posameznih delov brez poprejšnjega dovoljenja založnika ni dovo- ljeno. Dot'il'\lo' NlcT"Ll'" 1__ III IL _, II~ _ RAČUNANJE VREDNOSTI NEKATERIH MATEMATIČNIHFUNKCIJ Uvod Verjetno ste se že večkrat vprašali, kako kaikuiatorji ali pa tudi vaši domači ra- čunalniki računajo vrednosti različnih matematičnih funkcij. Če je funkcija dovolj "lepa", in take so prektično vse funkcije, ki jih srečamo v srednji šoli, potem lahko vrednost funkcije zapišemo v obliki neskončne vrste. Učeno bi lahko dejali, da funkcijo razvijemo v Taylorjevo vrsto . Tako na primer lahko zapišemo identiteto za funkcijo sin (x) : . x 3 x 5 x 7 sinlx] =x --+---+ ... 3! 5! 71 Neskončnega računa seveda ne moremo opraviti , toda če vzamemo za- dostno število členov vrste, lahko dobimo poljubno natančnen približek. Na prvi pogled smo tako problem računanja funkcij rešili , in to celo zelo splošno. Žal pa je v računalništvu, za razliko od matematike, ponavadi hitrost algoritma važnejši kriterij kot splošnost rešitve. Problem pri računanju vrednosti funkcij je s pomočjo Taylorjeve vrste je v tem, da moramo pri velikem argumentu x za dosego željene natančnosti upoštevati precej členov . To pa nasprotuje želji po učinkovitosti (hitrosti). Dobri algoritmi ponavadi zahtevajo globlje poznavanje problema, ki ga re- (y) 116 šujemo . Še preden si ogledamo, kako VeCII13 kalkulatorjev firme Hewlett- Packard računa kvadratni koren in trigonometrične funkcije, izpeljimo fo rmu- li, ki ju bomo potrebovali kasneje. Naj bo Tl poljubna točka v ravnini s koordinatama (Xl, vd. V polarnih koordinatah se koordinati točke Tl glasita: Vl = r. sin(\Pl) (1 ) kjer je r razdalja točke Tl od koordinatnega izhodišča in \Pl kot , ki ga poltrak iz izhodišča od Tl oklepa z osjo x. Zavrtimo točko Tl za kot \P2 v pozitivn i smeri okrog izhodišča v točko T2 (slika 1). Vprašajmo po koordinatah točke T2 • Najprej je: V2 = r . sin(\Pl + \P2) Če uporabimo za sin in cos adicijska teorema in upoštevamo ( t) . dobimo formuli : (2) Povejmo še, da kaikuiatorji HP uporabljajo posebne mikroprocesorje, pri katerih množenje z lan ne pomeni nič drugega kot premik decimalne pike. Prednost desetiške aritmetike je v tem, da ni potrebno pretvarjanje med dvo- jiškim in desetiškim sistemom. Če bi bila aritmetika dvojiška , bi računali dvo - jiške cifre, ki b i jih potem pomnožili z ustrezn imi potencami 2n • Kvadratni koren Kako bi izračunali kvadratni koren, če bi b ili brez kaikuiatorja, v usmerjeni šoli pa vas tega niso naučili? Izberimo neki približek, ki se nam zdi primeren, in prever imo , koliko sreče smo imeli , tako da izračunamo kvadrat približka. Če smo dob ili premalo , smo izbrali premajhen približek, če smo dob ili preveč, je bil približek prevelik , sicer pa smo imeli srečo in zadeli rešitev. Ker je verjet- nost zadetka majhna, postopek nadaljujmo tako , da približek izboljšamo, To delamo toliko časa, da smo z rezultatom zadovoljni. Zapišimo opisani postopek še v algoritmičnem jeziku. Vzemi začetni približek za koren a 117 ponavljaj izračunaj a2 R *"""x - a2 če je R dovolj majhen potem končaj sicer popravi a glede na predznak R do tod Seveda ta postopek še ni dober. Kalkulatorju je treba povedati, za koliko naj popravi e, razen tega pa zahteva vsak korak kar precej časa: izračun a2 in ostanka R. Algoritem torej spremenimo tako, da računanje a2 in R opravimo kar najlažje, ko spremenimo približek a. To naredimo tako, da postopoma računamo približek: najprej določimo pravilne tisočice, nato stotice, desetice, enice, .... Zdaj spoznamo v njem starega znanca (ali pa tudi ne). Tako vendar računamo kvadratni koren "peš" s papirjem in svinčnikom. Vpeljimo nekaj oznak, ki nam bodo olajšale nadaljnje delo: x število, katerega koren računamo a približek za koren b naslednja cifra X 1/2, ki jo iščemo j eksponent potence števila 10 Ra X - a2 tekoči ostanek aj novi e, ko dodamo b na ustrezno mesto: aj = a +b . 10i in Rb naj bo enak a/ - a2 Število x V2 bomo poiskali tako , da se bomo pravi vrednosti bližali od spodaj. Torej mora v vsakem trenutku veljati pogoj: in s tem Ra~O Ostanek moramo karseda zmanjšati, kljub temu pa še ostati pod X 1/2• Zato mora biti b največje tako število, a bo veljalo: 118 Če upoštevamo definicijo Rb Rb = (a + (b . 1Ql'))2 - a 2 = 2a b . 10i + (b . 10i)2 je b t isto največje število . da velja: 2ab . 10 i + (b . 10 i)2 .;;;; Ra Ko najdemo največje naravno število b, ki tej neenačbi zadošča, iz računa- mo nov približek a +- a+b .10 i Ra +- Ra - Rb j +- j - 1 Oglejmo si primer: x =54756, j = 1, a =200, Ra = 14756 . Zaporedoma računamo nove ostanke : b o O 1 4100 2 8400 3 12900 4 17600 14756 10656 6356 1856 -2844 Torej moramo za b vzeti 3, ker je b = 4 že preveč. Novi približek je tedaj a =200 + 3.101 =230, novi ostanek Ra = 1856 inj =O. Postopek nato ponovimo tolikokrat. da dosežemo željeno natančnost. Algoritem pa lahko tudi izboljšamo. Prva taka izboljšava je način , kako izrazimo (b . 10 i )2 . Pri tem upcštevamo .da je cšvscta prvih b lihih števil. Od tod: prištejmo 2a b . 10 i in dobimo 2ab. lO i + (b .lO i)2 = ~ (2a . ro/ + (2;-1) . 10 2i ) Leva stran je novi ostanek Rb . Namesto da gledamo Ra in Rb • opazujmo 5Ra in 5Rb • saj če velja Rb .;;;; Ra. velja tudi 5Rb ';;;; 5Ra. Ostanek Rb ima obliko 5Rb = ~ 10a . lO i + (10; - 5). 102i 119 postopek pa ostane nespremenjen, le da iščemo tisto največje naravno število b, za katerega velja 5R b .;;; Ra. novi ostanek pa je 5Ra =5Ra - 5R b . Trans- formacija je na prvi pogled nesmiselna. Če pa si ogledamo zaporedne Rb in upoštevamo, da zaradi BeD mikroprocesorja množenje z 10 pomeni le rotacij o v desno, dobi transformacija smisel. b=1 Rb=10a.10 j+05.10 2j b = 2 Rb = 10a. 10 j + 15 . 102j b = 3 Rb = 10a . 10 j + 25 . 102j Vidimo, da imajo ostanki Rb obliko 10a . lO j + ib - 1)15. 10 2j• kjer oznaka 15 pomeni dopisovanje petice. Poglejmo tako popravljeni postopek na našem primeru (x = 54756,j = 1, a = 200, 5Ra = 73780). Zaporedoma računamo nove ostanke: b 1 2 3 4 10a. 10 j + (b-1)15 .102j 20500 21500 22500 23500 53280 31780 9280 (novi 5Ra ) -14220 prekoračitev in naslednji korak lj = O. a = 230, 5Ra = 9280): b 1 2 3 4 5 10a.lOj + (b-1)\5.10 2j 2305 2315 2325 2335 2345 6975 4660 2335 O -2345 prekoračitev Če pogledamo obe tabeli, opazimo, da vrednosti v srednjem stolpcu ni treba računati, saj jih dobimo enostavno s spajanjem naračunanih cifer približ- ka, petice in še 2j ničel. Ničlo potem nadomeščamo zaporedoma s ciframi 1. 2, .... tako da je edino delo, ki ga še mora opraviti računalnik, računanje razlike 5R a - 5R b . Pri ogledu algoritma smo vzeli kot primer veliko število. Dejansko pa je v kalkulatorjih tipa HP štev ilo predstavljeno v eksponentnem zapisu z mantiso Min eksponentom exp: x =M. lOexp 120 kjer je 1 .;;;; M < 10. Poleg mantise je treba koreniti tudi potenco 1Oexp . S so- dim eksponentom ni t ežav: Če je eksponent l ih , pa ga zmanjšamo za 1. S tem smo seveda poveča li mantiso, tako da je v mejah 1 .;;;; M < 100, a po korenjenju bo mantisa M v pravih mejah. Ves postopek bi bi l torej takšen: 1. izračunaj eksponent odgovora 2. pomnoži mantiso s 5, dobiš 5Ra 3. z začetnim približkom a = O uporabi opisano metodo in po išči 12 cifer odgovora 4. zaokroži mantiso, dodaj eksponent 5. izpiši rezultat Če ste postopek razumeli, ga sprogramirajte v nekem višjem programske m jeziku, recimo v basicu na vašem domačem računalniku. Trigonometrične funkcije: sin(lp), cos(lp), tg(lp), ctg(lp) Za izračun vseh štirih trigonometričnih funkcij bomo uporabili isti algoritem, kar vsekakor pomeni bistven prihranek ROM-a. V vsakem primeru najprej izra- . čunamo tg(lp) ali pa ctg(lp), nato pa po potrebi iz njiju sin(lp) in coste) s f ormu- lama: sin(lp) = -----=--- (3) ± ctg(lp) 1 + ctg2 (lp) (4) V primerih , ko dobimo kot lp v drugih kotnih enotah kot radianih, ga najprej pretvorimo v radiane. Ker so trigonometrične funkcije periodične s periodo 21T, reduciramo poljuben kot na intervalu [O, 21T). To lahko naredimo tako , da odštevamo 21T od kota , dokler ne pridemo v ustrezn i interval. Vendar bi bil tak postopek za velike kote sila neučinkovit. Zato kot najprej zapišimo v eksponentnem zapisu in nato od njega odštevajmo velike večkratnike kota 21T, dokler ne dobimo negativnega kota. Nato ta večkratnik enkrat prištejemo in postopek ponovimo pri manjšem večkratniku kota 21T. Zapišimo ta postopek algoritmično: 121 kot 'Pzapiši v obliki 'P ==81.8283 .... 10n za k = n, n - 1, ..., Oponovi y~2.1T.1ok ponavljaj 'P~ 'P - Y dokler ni 'P< O 'P~'P+Y Na ta način velike kote zelo hitro spravimo v ustrezni integral. Sam pre- misli, kaj moramo opraviti v gornjem postopku, če imamo za podatek negati- ven kot. S tem algoritmom smo spravili kot na interval [O, 21T) in od tu dalje nas bodo zanimali le še ustrezni koti. Glavna ideja, ki nas pripelje do končnega algoritma, je naslednja. Če poznamo poljubno točko na poltraku iz izhodišča, ki oklepa z osjo x kot 'P, znamo tg('P) in ctg('P) izračunati s srednješolskima formulama: )1' tg('P) =- x in x ctg('P) =- Y (5) in odtod s formulama (1) in (2) tudi sin('P) in coste). Seveda pa še ne vemo , ka- ko naj pridelamo ustrezno točko (x, y) na poltraku, ki oklepa z osjo x kot 'P. Postopajmo takole: začetno točko, ki je dovolj blizu osi x, zavrtimo nekajkrat v pozitivni smeri za ustrezne kote 'P, dokler ne dobimo zaželjene točke . Najprej zapišimo formuli (2). ki zavrtita točko (Xl, Yd v ravnini za kot 'P2 : X2 = Xl . COS('P2) - YI. sin('P2) Y2 = YI . COS('P2) + Xl. sin('P2) Če obe enačbi (6) delimo s COS('P2), dobimo identiteto: X2/COS('P2)=XI-Yl.tg('P2) YdCOS('P2) == Y l + Xl· tg('P2) (6) Označimo prvo desno stran v (7) z X2', drugo desno stran pa z Y2', Če enačbi (7) zapišemo v teh novih oznakah, dobimo razmerje: (8) Ker je Y2/X2 = tg('Pl + 'P2), lahko tangens kota 'PI + 'P2 dobimo tako, da izra- čunamo X2' in Y2', Torej lahko izračunamo tangens za 'P2 večjega kota, če le poznamo tg('P2). X2' in Y2' To in pa seveda dejstvo, da postopek lahko večkrat po novimo, uporabimo v našem algoritmu. Za izračun tg('Pl + '(2) v (8) mora- mo izračunati X2' in Y2', ki ju dobimo v formuli (7). 122 Ker smo pri izbiri kota 1{)2 še svobodni. ga izberimo tako. da bo izračun v (7) čim preprostejši. Kot smo že omenili. uporabljajo HP kaikuiatorji posebne procesorje. pri katerih predstavlja množenje s potencami 10 le premik decimal- ne pike. Zato zahtevaj mo. da je tg(1{)2) oblike 1o-.Množenje v (7) namreč tedaj ni nič drugega kot premik decimaine pike za k mest. Kot 1{) torej zapišimo v obliki: Pri tem smo s tg-1 zapisal i inverzno funkcijo k funkcij i tangens. Vse kon, stante ao . a l •... so naravna števila. manjša ali enaka 10, tako da za vsako od njih potrebujemo en sam štiribitni zapis. Ustrezne približne kote zapišimo v radianih in stopinjah. tg-1 (1) tg-1 (0.1 l - 1 tg (0.01) tg-1 (0.001) 0.785398163 0.099668652 0.009999667 0.001000000 tg'l (1) tg'l (tU l tg-1 (0.01) tg'l (0.001) 45 5.710593137 0.572938698 0.057295604 Koti v radianih vsebujejo v svojih cifrah več pravilnosti. Vsi ti koti so seve, da stalno zapisani v ROM-u. večja pravilnost pa pomeni manjšo porabo prosto, ra. kar daje radianom prednost pred stopinjami. Zato smo tudi takoj na za- četku pretvorili kot v radiane. Razcep opravimo s preprostim algorimnom. Od kota odštevamo ustrezni kot. dokler kot ni negativen . potem pa mu še enkrat prištejemo isti kot. Ves postopek nato ponovimo na manjšem kotu . Zapišimo ta postopek zopet v obliki algoritma: odštevaj: za i = O. 1, 2. ... ponovi aj +- O ponavljaj 1{) +-1{) - tg-1 (10-j) aj +-aj + 1 če je 1{) < Opotem 1{) +-1{)+ tg'l (1 o-jI aj +-aj - 1 zapusti odštevaj do tod do tod 123 Kako majhen naj bo ostanek r pr i razcepu (9). je odv isno od natančnosti aritmetike. V večini kalkulatorjev HP se napravi razcep do kota tg-1 (0.0001) . Tedaj je namreč ostanek že tako majhen. da ne vpliva niti na zadnjo decimaiko v končnem izračunu. Da pa lahko začnemo celoten postopek. moramo določiti še začetno točko. Ker je ostanek r zelo majhen. je. če vzamemo x = 1, tg(r) približno enak rln za začetno točko vzamemo točko (1,r). Napišimo še algor i- tem, s katerim izračunamo končno točko na poltraku, ki oklepa z osjo x kot !(Jo x+-1 y +-r za k =O, 1...., m ponovi zai= 1.2• ...•ak ponovi x' +- x - y . lO-k y' +- Y + x . 10-k x -