PA081: Programování numerických výpočtů 2. Numerická stabilita (nejen) elementárních výpočtů Kvadratické rovnice Elementární funkce Algebraické úpravy Aleš Křenek Algebra kvaternionů jaro 2012 Kvadratické rovnice rovnice tvaru ► školní vzorec ax2 + bx + c = 0 -b ± yjb2 - 4ac 2a Kvadratické rovnice Elementárni funkce Algebraické úpravy Algebra kvaternionu Shrnuti Kvadratické rovnice rovnice tvaru školní vzorec ax + bx + c = 0 x+ -b ± -Jb2 - Aac 2a je-li b2 » ac, počítáme x+ z rozdílu dvou velmi blízkých čísel ► problém katastrofálního zrušení „»" neznamená až příliš velký rozdíl koeficientů ► typ f 1 oat - přesnost na 7-8 dekadických číslic ► při a = 1 tedy stačí b ~ 103c Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionů 2/35 Kvadratické rovnice konkrétně pro b = 2000, c = 1 ► „přesné" hodnoty ► VĎ= 1999.9989999997499998749999218749453... ► x+ = -.00050000012500006250003906252734377... ► x- = -1999.9994999998749999374999609374... ► pro f 1 oat ► VĎ= 1999.999 ► x+ = -0.00048828125, chyba 2.5% ► X- = -1999.9995, chyba odpovídá přesnosti typu Kvadratické rovnice ► numericky stabilní řešení pro x+ ► vlastnosti kořenů rovnice: axz + bx + c = (x - x+)(x - X-) tedy c = x+X- a lze počítat x+ = c/x_; ► pro předchozí příklad dostáváme ► x+ = -0.00050000014, chyba odpovídá přesnosti typu ► analogicky pro b < 0 je nepřesné x~, vypočte se symetricky Kvadratické rovnice ► numericky stabilní řešení pro x+ ► vlastnosti kořenů rovnice: axz + bx + c = (x - x+)(x - X-) tedy c = x+X- a lze počítat x+ = c/x_; ► pro předchozí příklad dostáváme ► x+ = -0.00050000014, chyba odpovídá přesnosti typu ► analogicky pro b < 0 je nepřesné x~, vypočte se symetricky ► i v triviálním výpočtu může vzniknout problém ► řešení může být docela jednoduché Elementární funkce goniometrické, hyperbolické, exponenciální, logaritmické, odmocniny, ... Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionú Elementární funkce goniometrické, hyperbolické, exponenciální, logaritmické, odmocniny, ... pohled matematika ► zaběhaný intuitivní aparát ► příjemné analytické vlastnosti (spojitost, derivace, ...) ► většina matematických modelů se bez nich neobejde Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionů Elementární funkce goniometrické, hyperbolické, exponenciální, logaritmické, odmocniny, ... pohled matematika ► zaběhaný intuitivní aparát ► příjemné analytické vlastnosti (spojitost, derivace, ...) ► většina matematických modelů se bez nich neobejde pohled programátora ► samy o sobě numericky stabilní ► optimalizované implementace v knihovnách ► problematické chování v okrajových případech vede na numerickou nestabilitu ► netriviální výpočetní náročnost Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionů ■o0.o 5/35 Elementární funkce Co je špatně? celý aparát vede na iracionální čísla (tt, e, VŽ,...) ► ve většině případů zdroj nepřesnosti ex ► tendence k přetečení - e88 = 1.65 x 1038 tanx ► přetečení pro x — j (a perioda) roste nade všechny meze sin x a cos x ► pro x — j resp. x — 0 vedou na špatně podmíněné rovnice ► malá změna vstupu vede k velké změně výstupu ► např. sinx = t pro t — 1 Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionů ■o0.o 6/35 Co s tím? ► uvědomit si možné problematické chování ► potřebuji skutečně tyto funkce pro svůj výpočet? ► mnoho problémů má i jiné řešení ► provést transformace, kterými se vyhneme numerické nestabilitě ► podobné (a + b)2 = a2 + 2ab + b2 ► podle kontextu zvolit vhodnou implementaci Algebraické úpravy výraz VI + x - -Ji - x pro malá x odčítání velmi blízkých čísel transformace X X (VI + x - VI - a:)(V1 + x + VI - a:) VI + x + VI - ^ (1 +x) - (1 -x) Vl + x + Vl - a: 2x VTTx + VI - x sčítání velmi blízkých čísel - dostatečně přesné Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionú Algebraické úpravy výraz ln -Jx + T - ln -fx pro velká x rozdíl blízkých čísel transformace Kvadratické rovnice Elementární funkce Algebraické úpravy x 1 - ln -Jx ln(x lnx2 1 (ln(x + 1) - lnx) 1, x + -ln- 2 x 1 ^ln(l + -) 2 x možná ztráta přesnosti, ale i tak lepší Algebra kvaternionů ■o0.o 9/35 Mocninné řady ► Taylorův rozvoj základních funkcí sin x = x cos x ln(l + x) = x X X X Tj + X2 X3 ~3!~ + ^! "' X2 2! f 4! "■ X2 x3 2 X x2 1 3x3 2 ~ 2-4 + 2 ■4-6 X 3x2 3 ■ 5x3 2 + 2 ■ 4 2 ■ 4 ■ 6 + 4 □ < Ö1 ► < i » Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionů 10/35 Mocninné řady Demonstrační příklad vyraz 1 - cosx X~2 se pro x —■ 0 blíží \ pro malá x odečtení dvou blízkých čísel konkrétně pro x = 0.0005 ve float: cosx = 0.99999988 1 - cosx = 1.1920928 x 10~7 (falešná přesnost) 1 ~C°SX = 0.47683709 xl správný výsledek je 0.49999999 (na 8 míst) Mocninné řady Demonstrační příklad náhrada cosx Taylorovým rozvojem X2 X4 x6 cosx = 1---1---- 2! 4! 6! pro výpočet na 8 cifer stačí do řádu x4 ? 4 1 1 1 i * * 1 ? 1-cosx 1 - 1 + "2--i2 1 x x~2 = x2 = 2 ~ 12 ► výsledek výpočtu pro x = 0.0005 je 0.49999997 ► podstatně lepší přesnost ► řádově méně aritmetických operací Mocninné řady Nedostatky ► rychlá konvergence pro malá x, pomalá pro větší ► přímý důsledek Taylorovy věty ► vyplatí se použít vlastní vzorec pro malé x ► resp. blízko problematického bodu numerické stability ► v ostatních případech zůstat u knihovní funkce ► jak poznáme, co je „malé x"? Kdy použít Mocninné řady ► známe požadovanou přesnost ► bezpečné je použít přesnost daného typu ► pro f 1 oat a výsledky - 1 je to 10"7 ► při méně přesných vstupních datech méně striktní ► nerovnost x < c musí zajistit, že největší zanedbaný člen řady je menší než požadovaná přesnost Kvadratické rovnice Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Mocninné řady Kdy použít známe požadovanou přesnost ► bezpečné je použít přesnost daného typu ► pro f 1 oat a výsledky ~ 1 je to 10~7 ► při méně přesných vstupních datech méně striktní nerovnost x < c musí zajistit, že největší zanedbaný člen řady je menší než požadovaná přesnost konkrétně v předchozím příkladu tedy x< 0.09212 zkontrolujeme, zdali dává smysl pro problematickou funkci cosO.09212 = 0.99999871 je ještě v toleranci □ ► Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionů 14/35 Mocninné řady Kdy použít výsledný kód pro 1 - cos x x2 Kvadratické rovnice Elementární funkce Algebraické úpravy float x,y; if (x < 0.09212) y = 0.5 + x*x/12.0; el se y = (l-cos(x))/(x*x); Algebra kvaternionú 15/35 Odmocnina ve většině používána k normalizaci vektorů silové působení, např. dva bodové náboje a a b ► velikost síly F = k, lalb F = F a-b kde r = ||a - b|| ► silový vektor snadný výpočet r2 = X(aŕ-fcŕ)2 r už vyžaduje odmocninu Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionů 16/35 Odmocnina Masivní použití počítačová grafika - osvětlení plochy ► zpravidla interpolací počítáme normály k zakřivené ploše ► vektory mají správný směr ale nejsou jednotkové ► pro správný výpočet osvětlení potřebná normalizace ke všem těmto výpočtům nepotřebujeme striktně -Jx -j= se hodí mnohem více ► násobení je lepší než dělení navíc zpravidla stačí velmi hrubá aproximace ► citlivý je hlavně směr vektoru, ten zůstává přímo Taylorův rozvoj -yj= iterační metody (Newtonova atd.) Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionú □ s 4 : ■oq.o 17/35 Odmocnina Fast inverse square root ► hra Quake III na konci 90. let ► převratně realistická grafika ► vděčila velmi rychlému výpočtu *■ kód zřejmě pochází z grafických knihoven SGI Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionů 4 □ ► 1); y = u.x * (1.5f - xhalf * u.x * u.x); return y; Kvadratické rovnice Elementárni funkce Algebraické úpravy Algebra kvaternionú Odmocnina Fast inverse square root ► funguje na základě IEEE reprezentace čísla x = (1 + mx) x 2ĚX III 1 0 0 1 1 1 1 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ■ 0.15625 31 30 23 22 (bit index) přitom mx = Mx/223 a ex = základ výpočtu 127 Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionú y \og2y log2(l + my x -log2x 1 log2(l + ra„ 19/35 Odmocnina Fast inverse square root ► pro m, G [0,1) si lze dovolit aproximaci log2(l + m) = m + a ► dosazením a úpravami dostaneme Ey X 223 + My R 1 (Exx2 23 tj. červený řádek v kódu s empiricky stanovenou konstantou R více viz http://en.wi ki pedi a.org/wi ki/Fast_ i nverse_square_root Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionú 20/35 Odmocnina Fast inverse square root poslední krok je jedna iterace Newtonovy metody pro vstup x hledáme takové y aby y = ^= tj. hledáme kořen funkce f(y) = -^-x Newtonovou metodou f(yn y-n+i = y n f (yn yn 1 x ^2 Vn yn yn(i což už odpovídá poslednímu výrazu v kódu Kvadratické rovnice Elementárni funkce Algebraické úpravy Algebra kvaternionú □ bude potřeba odmocnina Odmocnina Lze se jí i vyhnout výpočet velikosti síly F - — - -ľšA + — dr r13 r7 liché exponenty => bude potřeba odmocnina zajímá mě většinou silový vektor F = pa~b = (a - b) (— - — tedy vystačím s r2 Odmocnina Lze se jí i vyhnout ► implementace float n,r2,r4,r8,rl4,F[3]; i nt i ; foř (n=0,i=0; i<3; i++) { F[i] = a[i] - b[i]; n += F[i]*F[i]; } r2 = 1.0/n; r4 = r2*r2; r8 = r4*r4; rl4 = r8*r4*r2; n = 12*A*rl4 - 6*B*r8; f oř (i=0; i<3; i++) F [i] * Parametrizace rotace Aneb jak se vyhnout goniometrickým funkcím ve 2D jeden úhel

— qxq ortogonální lineární transformace Kvaterniony ► více ekvivalentních definic ► nejjednodušší M = K4 s kanonickou bází (1, i, j, k) a násobením i2 = f = kz = ijk = -1 ► chápeme ťcH jako ryze imaginární kvaterniony ► potom pro \q\ = 1 je zobrazení x >— qxq ortogonální lineární transformace ► nakrytí 2:1 (q a -q reprezentují stejnou transformaci) ► skládání rotací je násobení ► inverze je q ► korekce degenerace normalizací Interpolace rotací ► použití pro animace ve 3D ► lineární interpolace nestačí lerp(q0,qi,t) = (1 - t)q0 + tqi PA081: Programování numerických výpočtů A. Křenek I Kvadratické rovnice Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionú Interpolace rotací ► použití pro animace ve 3D ► lineární interpolace nestačí lerp(q0,qi,t) = (1 - t)q0 + tqi ► sférická lineární interpolace (SLERP) sin(l - t)Q sintQ slerp(q0, qi, t) = -—--q0 + . _ qi sinQ sinQ Pí f VI Kvadratické rovnice Elementární funkce Algebraické úpravy Algebra kvaternionú 30/35 Interpolace rotací ► použití pro animace ve 3D ► lineární interpolace nestačí lerp(q0,qi,t) = (1 - t)q0 + tqi ► sférická lineární interpolace (SLERP) slerp(q0,qi,t) sin(l - ř)Q sinQ 1o + sintQ sinQ pí ► v kvaternionech lze vyjádřit slerp(q0,qi,t) = qoiq^qi)1 Superpozice množiny bodů Formulace problému ► množiny odpovídajících si bodů {xí} a {yŕ} v K3 ► hledáme transformaci Q (rotaci + posunutí) tak, aby XllQxŕ-yŕlI2 bylo minimální ► posunutí - stačí nezávisle posunout těžiště {xí} do těžiště {yr} ► zbývá najít rotaci Superpozice množiny bodů Použití kvaternionů formulace pomocí kvaternionů Xllqxŕíí-yŕ||2 = X(llxŕ|l2 + 2 ^Re(yŕqxŕ