Kako se zaista računaju sin(), cos(), log(), exp() itd, ukratko

Kod digitalne obrade signala u realnom vremenu kritično je važna brzina izvršavanja programskih funkcija koje računaju koren, sinus, kosinus, logaritam itd.

A vrednosti mnogih važnih matematičkih funkcija u principu nije moguće izračunati – u konačnom vremenu – uz pomoć četiri osnovne računske operacije, korenovanja i stepenovanja.

Ako se neko oduvek pitao kako se onda one uopšte izračunavaju u/na silicijumu…

To računanje načelno je najbrže u namenskim modulima hardvera, koji na nivou arhitekture čipa primenjuju algoritme kao što je CORDIC i njegovi derivati. Ti algoritmi daju aproksimacije, tj. određuju približne vrednosti pomenutih matematičkih funkcija, ali takve da je greška po pravilu zanemariva. CORDIC koristi binarnu pretragu i iteracije zasnovane na suptilnoj (mada ne toliko složenoj) matematici, koja je omogućila i to da usputne uzastopne operacije budu svedene na pomeranje bitova i sabiranje, što se lako primenjuje i brzo izvršava na nivou hardvera.

Mnoge hardverske arhitekture nemaju ugrađenu tu mogućnost (mada vrlo napredni mikrokontroleri imaju makar hardversku jedinicu za računanje s pokretnim zarezom, koja ubrzava računanje generalno, ali to je druga tema).

Ako posmatramo samo mogućnosti hardvera, najskromniji tipovi osmobitnih mikrokontrolera umeju samo da sabiraju i oduzimaju cele brojeve, a većina mikrokontrolera nudi još i hardversko množenje. (U ovu drugu grupu spada i ATmega328p u Arduinu Nano, Uno itd.)

Zato se prilagođeni algoritmi za aproksimaciju vrednosti raznih matematičkh funkcija primenjuju u softveru, pre svega u kompajlerima, tj. unutar funkcija programskih jezika. Iza sin(), cos(), log() itd stoje algoritmi kojima se do dovoljno tačnog rezultata dolazi uz upotrebu samo osnovnih računskih operacija, i to – vodeći mnogo računa o efikasnosti, tako da se za dobijanje rezultata što veće tačnosti potroši što manje vremena.

Ti algoritmi često su zasnovani na polinomima Čebiševa, a odabrani su i pažljivo optimizovani prema mogućnostima konkretne hardverske arhitekture. Tu vidimo matematičke izraze koji na prvi pogled nemaju nikakve veze sa – na primer – sinusom, ali nekom matematičkom “magijom” dobije se rezultat koji je dovoljno blizu vrednosti sinusa datog ugla. Teorijski, mogao bi se tako dobiti i idealno tačan rezultat, ali bi za to bilo potrebno beskonačno vreme, a digitalno predstavljeni brojevi ionako imaju ograničenu preciznost – pa se bira kompromis.

Ovaj izraz, na primer, daje rezultat sin(x) s maksimalnom greškom od ±0,0000006 za x između -π/2 i π/2, a zahteva samo osnovne računske operacije. Polinom nižeg stepena (s manje članova) izvršavao bi se brže, ali bi greška bila veća, i obrnuto. Takođe ne treba zaboraviti da nesavršenost samog polinoma ovde nije jedini izvor greške. Svaka operacija s brojevima s pokretnim zarezom dodaće još i grešku zaokruživanja.
Apsolutna greška aproksimacije za gornji polinom. Talasast oblik grafika greške karakterističan je za polinome Čebiševa i nije posledica periodičnosti funkcije koja se aproksimira, ovde sin(x).
Radi poređenja: kod aproksimacije Tejlorovim polinomom greška bi bila vrlo mala oko jedne tačke i brzo bi rasla već na maloj udaljenosti od nje. Polinom Čebiševa (istog stepena) daje manju maksimalnu grešku za ceo dati raspon (interval).

Drugi primer, druga funkcija (polinom šestog stepena, maksimalna greška ±6*10-6 za -1 ≤ x ≤ 1):

ex ≈ 0.99999981 + 1.00002170x + 0.50000620x2 + 0.16649306x3 + 0.04163566x4 + 0.00868056x5 + 0.00143849x6

Aproksimacije imaju predvidivu i (recimo) prihvatljivu maksimalnu grešku samo unutar izvesnog raspona (intervala) vrednosti argumenta x. Pri izračunavanju, argument se prvo mora “uterati” u taj raspon, ako već nije u njemu. Na primer, ex, tj. exp(x), može biti prvo izraženo u obliku 2n * er, gde je n ceo broj, a r je unutar predviđenog raspona. Naravno, i taj račun već je deo funkcije u kompajleru.

Zelena linija je funkcija ex, a crvena je njena aproksimacija polinomom trećeg stepena, dakle, kraćim i bržim od prethodnog, uz veću maksimalnu grešku i unutar suženog raspona.

Osim ovakvih polinoma, za dobijanje “skoro tačnih” vrednosti raznih matematičkih funkcija mogu postojati i drugi izrazi i algoritmi, ponekad specifični “trikovi”, ponekad rekurzivni (kod generisanja sinusnog signala), a upotrebljivost svakog od njih zavisi od konkretne situacije. Pitanje se najčešće svodi na izbor između brzine i tačnosti. Ako, na primer, namena projekta zahteva brzinu, a može podneti grešku veću od one koju nudi ugrađena funkcija programskog jezika, ta funkcija može biti zamenjena kraćom i bržom aproksimacijom u programu.

Zelena linija: funkcija ex. Crvena: njena aproksimacija jednostavnim izrazom
((x+3)2 + 3) / ((x-3)2 + 3). (Vidi pojam Padeova aproksimacija.)
Za interval -0,5 ≤ x ≤ 0,5 greška je manja od 0,000075.

Slika: hackaday.io

Načelno drukčiji put za dobijanje rezultata jeste upotreba tabele (lookup table) u kojoj je navedeno, na primer, 90 unapred izračunatih vrednosti sinusa za svih 90 celih stepeni jednog kvadranta. Simetrije omogućavaju da se odatle lako dobiju vrednosti sinusa i kosinusa i za sve ostale celobrojne vrednosti ugla. Ako je potrebno, za vrednosti između celih brojeva može se primeniti linearna interpolacija, što je brzo, mada će se opet dobiti (dovoljno) približna vrednost. Mana je to što tabela/tabele zauzimaju prostor u nekoj od memorija, ali svakako je brže povući vrednost iz niza nego računati. U sledećem primeru predstavljena je takva funkcija primenjiva na Arduinu, sa ili bez interpolacije, višestruko brža od ugrađene funkcije sin().
https://forum.arduino.cc/t/a-faster-sin-lookup-function/69173

Arduino za sin() koristi polinom s koeficijentima koji upućuju na Maklorenov (Tejlorov) niz. U izvornom kodu funkcije navedeni su koeficijenti do 11. stepena, što odgovara maksimalnoj grešci aproksimacije u redu veličine 10-8 za argument između -π/2 i π/2. Kosinus se računa preko iste rutine, samo se prethodno argument rotira za π/2. Inače, u Arduinu na 16 MHz izvršenje sin() traje tipično 110…120 μs, tan() tipično 130…170 μs. Izvorni kod u asembleru, i za ostale funkcije:
https://github.com/avrdudes/avr-libc/tree/main/libm/fplib

Dalje, razno:
https://www.embeddedrelated.com/showarticle/152.php
https://scholarworks.umt.edu/cgi/viewcontent.cgi?article=1313&context=tme
https://wiki.logre.eu/index.php/Trigonom%C3%A9trie_en_virgule_fixe
https://www.youtube.com/watch?v=4SsspqSdwaA

https://namoseley.wordpress.com/2012/09/21/a-better-sinx-approximation/

https://bmtechjournal.wordpress.com/2020/05/27/super-fast-quadratic-sinusoid-approximation/

http://www.ganssle.com/item/approximations-c-code-exponentiation-log.htm

Leave a Reply