O čem to bude PA081: Programování numerických výpočtů Aleš Křenek Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy jaro 2015 Algebra kvaternionů O čem to bude zajímá nás chování skutečného světa ► problémy přírodovědné, technické, humanitní... popisujeme matematickými prostředky ► zejména pomocí reálných čísel ► umělý aparát, leckdy zcela neodpovídá skutečnosti ► za staletí celkem zvládnutý, vyučovaný od základní školy, obecně přijímaný Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionů O čem to bude ► zákonitosti zkoumaných systémů typicky vyjádřeny rovnicemi ► zajímavé systémy — složité rovnice ► reprezentativní příklad - Schrôdingerova rovnice ifr^-Y(x,í) =/TF(x,í) ot ► analyticky řešitelná jen pro triviální případy ► izolovaná částice, částice v potenciálové jámě, ..., atom vodíku ► i tak je řešení dost komplikované ► numerické řešení je jediný realisticky možný přístup O čem to bude ► numerická matematika ► řešení matematicky formulovaného problému aritmetickými prostředky ► zpravidla algoritmický postup - numerická metoda ► disciplina podstatně starší než počítače O čem to bude ► numerická matematika ► řešení matematicky formulovaného problému aritmetickými prostředky ► zpravidla algoritmický postup - numerická metoda ► disciplina podstatně starší než počítače ► metody jsou známé, naprogramujeme je a je hotovo ► ne tak docela (De)motivační příklad ► uměle zkonstruovaný, ilustruje podstatu problému ► pro dané n vypočítejte integrál f. Jo xnex ľdx Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionů (De)motivační příklad uměle zkonstruovaný, ilustruje podstatu problému pro dané n vypočítejte integrál f. Jo xnex 1dx očekávané vlastnosti ► E0 = [e*-i]l = l-\ ► pro 0 < x < 1 platí xn > 0 a 0 < e*'1 < 1, tedy En > 0 ► podobně En < xndx =-- a tedy lim En = 0 Jo n + 1 n-oo (De)motivační příklad pro numerický výpočet lze transformovat na rekurentní posloupnost integrací per-partes u = xn, dv = ex~ldx nxn xex 1dx 1 - nEn-i Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionů 4 □ ► 4 = ► < Problémy konečné reprezentace Násobení a dělení ► samo o sobě nezanáší významnou chybu ► zachovává existující chybu ► příklad: (a + b)2 pro a = 1.23, b = 0.0155, na 3 cifry ► přesný výsledek je 1.55127025 Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionů Problémy konečné reprezentace Násobení a dělení ► samo o sobě nezanáší významnou chybu ► zachovává existující chybu ► příklad: (a + b)2 pro a = 1.23, b = 0.0155, na 3 cifry ► přesný výsledek je 1.55127025 ► přímý výpočet na 3 platné cifry: a + b = 1.23 + 0.02 = 1.25 tedy {a + b)2 = 1.56, chyba 0.56% ► není to tak zlé, ale může být lepší Čísla v plovoucí řádové čárce Numerická stabilita Elementárni funkce Algebraické úpravy Algebra kvaternionú 15/67 Problémy konečné reprezentace Násobení a dělení po transformaci (a + b)2 = a2 + 2ab + b2: 1.51 + 2 x 0.0191 + 0.000240 = 1.51 + 0.0382 + 0.000240 = 1.55 chyba 0.082 %, tedy téměř o řád menší za cenu 6 aritmetických operací místo 2 bude 3x pomalejší Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú Problémy konečné reprezentace Násobení a dělení po transformaci (a + b)2 = a2 + 2ab + b2: 1.51 + 2 x 0.0191 + 0.000240 = 1.51 + 0.0382 + 0.000240 = 1.55 chyba 0.082 %, tedy téměř o řád menší za cenu 6 aritmetických operací místo 2 bude 3x pomalejší ... uvidíme Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú Měření rychlosti výpočtu Naivní přístup POSDÍ volání getti meofdayQ getti meofdayC&start,NULL) ; c = a + b; c *= c; getti meofdayC&stop,NULL); současné CPU 2-3 GHz 1 aritmetická operace ~ 1-10 ns nemáme tak přesné hodiny □ ô1 4 : Čísla v plovoucí řádové čárce Numerická stabilita Elementárni funkce Algebraické úpravy Algebra kvaternionú 17/67 Měření rychlosti výpočtu Zopakujeme v cyklu ► opakování 109x navýší čas na měřitelnou getti meofdayC&start,NULL) ; for (i=0; i<1000000000; i++) c = a + b; c *= c; } getti meofdayC&stop,NULL); Měření rychlosti výpočtu Zopakujeme v cyklu ► opakování 109x navýší čas na měřitelnou getti meofday(&start,NULL) ; for (i=0; i<1000000000; c = a + b; c *= c; } getti meofday(&stop,NULL); ► počítá podstatně rychleji ► trochu podezřelé Měření rychlosti výpočtu Zopakujeme v cyklu ► optimalizující kompilátor gcc -03 -f unrol 1-1 oops movss 44(%rsp), %xmm0 xorl %eax, %eax addss 40(%rsp), %xmm0 mul ss %xmm0, %xmm0 addl $8, %eax cmpl $1000000000, %eax j ne .L2 movss %xmm0, 36(%rsp) v těle cyklu se nic nepočítá optimalizaci obecně chceme ► použití registrů pro proměnné, max. využití FPU, eliminace zbytečného opakování je přilij ^ Měření rychlosti výpočtu Zopakujeme v cyklu ► brzdící kód ► uměle vložený do těla cyklu ► za běhu programu se nevyvolá - nebrzdí doopravdy ► zabrání příliš agresivní optimalizaci cyklu nikdy = (strlen(argv[0]) == 0) for (i=0; i<1000000000; { c = a + b; c *= c; if (nikdy) brzda(&nikdy,&a,&b,&c); Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú kompilátor musí předpokládat: ► proměnná ni kdy nemusí být 0 ► funkce brzdaQ má vliv na odkazované proměnné □ 4 S ► 4 = > 4 ^ * Měření rychlosti výpočtu Jak to dopadlo Intel Í7-2600 3.4 GHz vzorec čas Mflop/s Mcyklů/s (a + b)2 0.61 3278 1639 a1 + 2ab + bz 0.99 6060 1010 rozvinutý výpočet jen 1.6x pomalejší ► vnitřní paralelismus procesoru ► více FPU jednotek ► pipelining ► spekulativní výpočet větví programu ► stále dosahujeme jen 25 % FPU výkonu jednoho jádra CPU Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú Problémy konečné reprezentace Násobení a dělení Snažme se vzorce transformovat tak, aby se nejdříve násobilo/dělilo, pak teprve sčítalo/odčítalo. Je šance na přesnější výsledek, dopad na výkon nemusí být významný. Problémy konečné reprezentace Přetečení a podtečení ► násobení dvou velkých čísel může vést k nereprezentovatelnému exponentu 1.2e30 x 1.2e30 = 1.44e61 ► typ float umí exponent [-37, 38] ► výsledek operace už je oo Problémy konečné reprezentace Přetečení a podtečení ► násobení dvou velkých čísel může vést k nereprezentovatelnému exponentu 1.2e30 x 1.2e30 = 1.44e61 ► typ float umí exponent [-37, 38] ► výsledek operace už je oo ► podobně násobení dvou malých čísel vede k podtečení ► podle nastavení aritmetiky je výsledek ±0 nebo denormalizované číslo O.OOOOOOmmm x 1(T37 ► další výpočty nepřesné a navíc citelně pomalejší Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú 23/67 Problémy konečné reprezentace Přetečení a podtečení Při násobení a dělení více čísel uspořádejme posloupnost operací, abychom nedělili velmi malé číslo velkým, velké malým, a nenásobili vzájemně velká nebo malá čísla. Problémy konečné reprezentace Porovnávání PA081: Programování numerických výpočtů A. Křenek H O čem to bude ► nevinný výpočet |číslav plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionú Shrnutí Problémy konečné reprezentace Porovnávání nevinný výpočet float a=1.505, bl=1.315, b2=1.695, c=a+a, d=bl+b2; printf("%f %f %s\n'\c,d,c == d ? "jsou stejná" : "jsou ruzna"); má překvapivý výstup 3.010000 3.010000 jsou ruzna Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionů Problémy konečné reprezentace Porovnávání ► nevinný výpočet float a=1.505, bl=1.315, b2=1.695, c=a+a, d=bl+b2; printf("%f %f %s\n",c,d,c == d ? "jsou stejná" : "jsou ruzna"); má překvapivý výstup 3.010000 3.010000 jsou ruzna přesnější formát výpisu "%15 .12f": 3.009999990463 3.010000228882 jsou ruzna Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú 25/67 Problémy konečné reprezentace Porovnávání ► operátor == na reálných číslech nemá valný smysl ► téměř vždy je zasažen chybou předchozího výpočtu ► místo něj test na dostatečnou blízkost fabs(c-d) < EPSILON ► stanovení toleranční konstanty zpravidla empiricky ► různá pro hodnoty le-15 a le30 ► silné záleží na dané úloze ► ovlivněna i každým konkrétním výpočtem Problémy konečné reprezentace Porovnávání ► porovnávání < a > trpí stejným problémem ► záludnější podoba a komplikovanější řešení ► chci porovnat přesné hodnoty x a y ► znám jen přibližná (spočtená) x = x ± ex, ý = y y Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy y Algebra kvaternionú Problémy konečné reprezentace Porovnávání ► porovnávání < a > trpí stejným problémem ► záludnější podoba a komplikovanější řešení ► chci porovnat přesné hodnoty x a y ► znám jen přibližná (spočtená) x = x ± ex, ý = y ± ey: i-•-1 x y y opatrný („miserly") přístup ► intervaly se i částečně překrývají => tvrdíme x = y ► jinak porovnáme x ^ ý Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú Problémy konečné reprezentace Porovnávání ► porovnávání < a > trpí stejným problémem ► záludnější podoba a komplikovanější řešení ► chci porovnat přesné hodnoty x a y ► znám jen přibližná (spočtená) x = x ± ex, ý = y ± ey: i-•-1 x y y opatrný („miserly") přístup ► intervaly se i částečně překrývají => tvrdíme x = y ► jinak porovnáme x ^ y hladový („eager") přístup ► chceme vědět, že mohlo nastat x > y porovnávame x + ex > y - e y- Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú 27/67 Problémy konečné reprezentace Porovnávání Na rovnost porovnávejme vždy s tolerancí. U porovnání na nerovnost si uvědomme, co chceme vědět. Příklad Kvadratická rovnice rovnice tvaru školní vzorec ax + bx + c = 0 x+ = -b ± \/b2 - 4ac 2a Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra k\ atemionů 4 □ ► 4 f5> ► < í ► Příklad Kvadratická rovnice ► rovnice tvaru ► školní vzorec ax + bx + c = 0 -b ± \lb2 - 4ac 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 Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú 29/67 Příklad 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 Příklad Kvadratické rovnice ► bezproblémové ř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 Příklad Kvadratické rovnice ► bezproblémové ř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é Numerická stabilita (Pseudo)definice ► „metoda počítá sice špatně, ale jenom trochu špatně" ► žádoucí a přitom realistická vlastnost všech numerických algoritmů Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú 4 □ ► 4 9 k 4 ■= * 4 -E Numerická stabilita (Pseudo)definice ► „metoda počítá sice špatně, ale jenom trochu špatně" ► žádoucí a přitom realistická vlastnost všech numerických algoritmů ► pro vstup x hledáme řešení y = f (x) ► ve skutečnosti ► na aproximaci vstupu x = x + ex *■ nepřesnou numerickou metodou / ► dostaneme výsledek y = y + ey Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú 4 □ ► -7 zkontrolujeme chování výpočtu v c (v double): cosc = 0.99576087237414645243 (1 - cosc)/c2 = 0.499646 58945642923198 1/2 - c2/24 = 0.49964648949583333334 c4/6! = 0.000000 09997574124493 Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú Mocninné řady Kdy použít výsledný kód pro float x,y; 1 - cos x X2 if (x < 0.09211) y = 0.5 + x*x/24.0; el se y = (l-cos(x))/(x*x); Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionů 4 □ ► < fl> ► < -š ► < -Š ► 46/67 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 kde r = ||a - b|| silový vektor r ^ a-b F = F ■ - snadný výpočet r2 = X(aŕ-fcŕ)2 r už vyžaduje odmocninu 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.) □ s 4 : ■OQ.O Čísla v plovoucí řádové čárce Numerická stabilita Elementární íunkce Algebraické úpravy Algebra kvaternionú 48/67 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 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 Odmocnina Fast inverse square root ► funguje na základě IEEE reprezentace čísla x = (1 + mx) x 2ĚX _significant! (23 bits)_ sign exponent (8 bits) 11- ■ 0.15625 31 30 23 22 přitom mx = Mx/2Zi a ex základ výpočtu (bit index) Ei 127 Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy y \og2y x log2(l + ra3 -log2x 1 log2(l + mx Algebra kvaternionú 50/67 Odmocnina Fast inverse square root ► pro m, G [0,1) si lze dovolit aproximaci log2(l + m) = m + a ► dosazením a úpravami dostaneme 1 Ey x 2 5 + My = R --(Ex x 2Zá + Mx) 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 Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú 51/67 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 iy) = ^2 - 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 Čísla v plovoucí řádové čárce Numerická stabilita Elementární 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; { 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; F [i] * Odmocnina Lze se jí i vyhnout ► implementace float n,r2,r4,r8,r 14,F[4]; i nt i ; foř (n=0,i=0; 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; for (i=0; i< ; F [i] *= n; Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú kvůli vektorizaci může mít smysl F [4] ■O0.O Parametrizace rotace Aneb jak se vyhnout goniometrickým fu ► ve 2D jeden úhel

— qxq rotace ve 3D Kvaterniony ► algebraický aparát na vektory ve 4D ► 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 rotace ve 3D ► nakrytí 2:1 (q a -q reprezentují stejnou transformaci) ► skládání rotací je násobení ► inverze je q ► korekce degenerace normalizací Kvaterniony součást mnoha knihoven ► libquat, Boost, Ogre, DirectX, ... ► optimalizované implementace typické funkce ► norma/normalizace ► komplement ► násobení ► transformace vektoru/vektorů ► převod z/na Eulerovy úhly převod z/na matice 3x3 v praxi zpravidla použijeme příhodnou knihovnu ► tyto operace už před námi někdo naprogramoval lépe ► musíme chápat principy, na kterých stojí Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú Interpolace rotací použití pro animace ve 3D lineární interpolace lerp(q0,qi,t) = (1 - t)q0 + tqx 1(1 - t)q0 + tqi\ nevhodné vlastnosti ► nelineární vztah t a úhlu rotace ► špatná kontrola úhlové rychlosti Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú Interpolace rotací sférická lineární interpolace (SLERP) sin(l - ř)Q sintQ slerp(qo,qi,t) = sinQ lo + —ň-fli kde cos Q = qo ■ 4i PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Interpolace rotací sférická lineární interpolace (SLERP) sin(l - t)Q sintQ slerp(q0, qi, t) = -—-q0 + . _ qi sin Q sin Q kde cosQ = qo ■ qi v kvaternionech lze vyjádřit slerp(q0,qi,t) = q^q^qxÝ velmi efektivní implementace s použitím SSE, AVX Interpolace rotací ► sférická lineární interpolace (SLERP) sin(l - t)Q slerp(q0,1i,t) kde cos Q = qo ■ 1i v kvaternionech lze vyjádřit sin Q -1o sinřQ sin Q li slerp(q0,ql3t) = qo(q01qi)t velmi efektivní implementace s použitím SSE, AVX pro většinu operací s rotacemi jsou kvaterniony nejlepší volba ► přítel Google ochotně poradí Čísla v plovoucí řádové čárce Numerická stabilita Elementární funkce Algebraické úpravy Algebra kvaternionú 66/67 Shrnutí matematické modely reality použití idealizovaného aparátu reálných čísel konečná reprezentace čísel je zdrojem problémů ► různé projevy pro různé operace ► i v triviálním výpočtu může dojít k numerickému problému ► míra závadnosti formalizována pojmem numerické stability elementární funkce ► samotná implementace numericky stabilní ► použití ve výrazech může vést na numerické problémy ► explicitní Taylorův rozvoj pro okrajové případy rychlost kritické části výpočtu ► kompromis rychlost vs. přesnost ► změření nemusí být triviální □ s 4 : ■OQ.O Čísla v plovoucí řádové čárce Numerická stabilita Elementární íunkce Algebraické úpravy Algebra kvaternionú 67/67