A. Křenek O čem to bude PA081: Programování numerických výpočtů Aleš Křenek jaro 2019 Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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ý Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 Íft|-Y(x,r) = HY(x,t) 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 PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí (De)motivační příklad ► uměle zkonstruovaný, ilustruje podstatu problému ► pro dané n vypočítejte integrál En — rl xnex ldx o PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí (De)motivační příklad ► uměle zkonstruovaný, ilustruje podstatu problému ► pro dané n vypočítejte integrál rl En — xnex ldx o ► očekávané vlastnosti ► E0 = [e^]l = l-l ► pro 0 < x < 1 platí xn > 0 a 0 < ex~x < 1, tedy En > 0 ► podobně En ^ xndx = o -- a tedy lim En = 0 U + 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 En= [ - *xnex 1 -a o nx n-l 0 dx = 1 nEn-i PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 6/72 (De)motivační příklad ► pro numerický výpočet lze transformovat na rekurentní posloupnost ► integrací per-partes u = xn, dv = ex~ldx En = xne -.1 _ 0 rl 0 nxn lex ldx = 1 nEn-i PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů (De)motivační příklad PA081: Programování numerických výpočtů i : 0.367879 A. Křenek 2 : 0.264241 3 : 0.207277 0 čem to bude 4 : 0.170893 Čísla v plovoucí 5 : 0.145534 řádové čárce 6. : 0.126796 Numerická stabilita 7: 0.112430 Zaokrouhlování 8: 0.100563 Elementární 9: 0.094933 funkce 10 11 12 13 14 15 16 17 18 19 0.050674 0.442581 -4.310974 57.042664 -797.597290 11964.958984 -191438.343750 3254452.750000 -58580148.000000 1113022848.000000 Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 7/72 (De)motivační příklad ► co je špatně? ► pracujeme s konečnou reprezentací reálného čísla ► společný problém ručního i strojového zpracování ► v počítači daleko plíživější podoba (nevidíme mezivýsledky) PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí (De)motivační příklad ► co je špatně? ► pracujeme s konečnou reprezentací reálného čísla ► společný problém ručního i strojového zpracování ► v počítači daleko plíživější podoba (nevidíme mezivýsledky) ► v proměnné E [n] není uloženo přesně En už na začátku zatíženo chybou, E [0] = Eq + e PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí (De)motivační příklad ► co je špatně? ► pracujeme s konečnou reprezentací reálného čísla ► společný problém ručního i strojového zpracování ► v počítači daleko plíživější podoba (nevidíme mezivýsledky) ► v proměnné E [n] není uloženo přesně En už na začátku zatíženo chybou, E [0] = Eo + e ► i při zcela přesném výpočtu: E [1] = 1 - E [0] = l-Eo-e = E1-e E[2] = 1 - 2E[1] = 1 - 2Ei + 2e = E2 + 2e E[3] = 1 - 3E[2] = 1 - 3E2 - 3(2e) = E3 - 3(2e) PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí E[n] = = En + {-l)nn\e O čem to tedy bude ► představení numerických metod pro řešení vybraných problémů ► pragmaticky, bez rozsáhlé teoretické analýzy, okrajových podmínek atd. ► formulace přiměřeně přesného a efektivního algoritmu ► matematicky korektní metoda nestačí ► pro různá použití jsou vhodné různé alternativy ► použití na smysluplném příkladu PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Čísla v plovoucí řádové čárce ► standard IEEE 754 ► vychází návrhu reprezentace čísel v koprocesoru Intel 8087 ► základní formát ± 1 .mmmmmmmmm... x 2±ee- ► znaménko mantisy (1 bit) ► mantisa l.mmmmmmmm ► binárně, absolutní hodnota ► číslo v intervalu [1,2) v dané přesnosti ► nekonečná množina pokryta konečným počtem hodnot ► základní zdroj nepřesnosti ► exponent ► binární číslo v rozsahu 1 až např. 254 ► znamená hodnoty exponentu -126 až+127 ► speciální význam hodnot 00... 00 a 11... 11 - kódování ±0, ±oo, ±NaN PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Čísla v plovoucí řádové čárce ► nejběžnější typy - velikosti v bitech a rozsah exponentu typ exp. mantis a celkem rozsah exp. Single (float) 8 23 32 127 Double 11 52 64 1023 Quad (long double) 15 112 128 16383 ► tj. např. největší číslo typu float je 2127 = 10 ► nej menší nenulové 2~126 = 10"37 ► relativní chyba je omezena 2~25 = 10"8 ► tedy počítáme na cca. 8 platných cifer Problémy konečné reprezentace Ztráta přesnosti sčítání ► pro zjednodušení v desítkové soustavě a na 3 platné cifry mantisy ► sečtěme 1.23e3 a l.OOeO ► nutné sjednocení exponentů: 1.23e3 + 0.001e3 = 1.231e3 = 1.23e3 ► v rámci dané přesnosti korektní výsledek PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Problémy konečné reprezentace Ztráta přesnosti sčítání ► pro zjednodušení v desítkové soustavě a na 3 platné cifry mantisy ► sečtěme 1.23e3 a l.OOeO ► nutné sjednocení exponentů: 1.23e3 + 0.001e3 = 1.231e3 = 1.23e3 ► v rámci dané přesnosti korektní výsledek ► k 1.23e3 desetkrát přičtěme l.OOeO ► opakované aplikování předchozího postupu dává opět 1.23e3 ► nemusí to být to, co jsme chtěli ► 1.23e3 + (l.OOeO + ■ ■ ■ + l.OOeO) = 1.24e3 PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Problémy konečné reprezentace Ztráta přesnosti sčítání ► pro zjednodušení v desítkové soustavě a na 3 platné cifry mantisy ► sečtěme 1.23e3 a l.OOeO ► nutné sjednocení exponentů: 1.23e3 + 0.001e3 = 1.231e3 = 1.23e3 ► v rámci dané přesnosti korektní výsledek ► k 1.23e3 desetkrát přičtěme l.OOeO ► opakované aplikování předchozího postupu dává opět 1.23e3 ► nemusí to být to, co jsme chtěli ► 1.23e3 + (l.OOeO + ■ ■ ■ + l.OOeO) = 1.24e3 ► operace v plovoucí řádové čárce nejsou asociativní ► ani tam, kde jsme na to z algebry zvyklí PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Problémy konečné reprezentace Katastrofální zrušení ► odečtení dvou vzájemně blízkých čísel vede k výrazné ztrátě přesnosti PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 13/72 Problémy konečné reprezentace Katastrofální zrušení ► odečtení dvou vzájemně blízkých čísel vede k výrazné ztrátě přesnosti ► rozdíl dvou měření s přesností na 3 platné cifry: 1.23 - 1.22 ► binární reprezentace: 1.0011101 a 1.0011100 ► zatížena chybami 0.3 % a 0.1 % Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Problémy konečné reprezentace Katastrofální zrušení ► odečtení dvou vzájemně blízkých čísel vede k výrazné ztrátě přesnosti ► rozdíl dvou měření s přesností na 3 platné cifry: 1.23 - 1.22 ► binární reprezentace: 1.0011101 a 1.0011100 ► zatížena chybami 0.3 % a 0.1 % ► odečtením dostaneme binárně 0.0000001, tj. 0.0078125 ► správný výsledek byl 0.01, relativní chyba 22 % ► navíc působí dojmem falešné přesnosti Problémy konečné reprezentace Katastrofální zrušení ► odečtení dvou vzájemně blízkých čísel vede k výrazné ztrátě přesnosti ► rozdíl dvou měření s přesností na 3 platné cifry: 1.23 - 1.22 ► binární reprezentace: 1.0011101 a 1.0011100 ► zatížena chybami 0.3 % a 0.1 % ► odečtením dostaneme binárně 0.0000001, tj. 0.0078125 ► správný výsledek byl 0.01, relativní chyba 22 % ► navíc působí dojmem falešné přesnosti Vyhněme se odčítání dvou blízkých čísel. Je-li to i tak nezbytné, počítejme s výsledkem zatíženým potenciálně velkou chybou. PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Problémy konečné reprezentace Robustní aritmetiky ► uvažujme přesnost na 6 cifer ► součet 10000.0 + 3.14159 + 2.71828 ► přesně 10005.85987, zaokrouhleně 10005.9, naivně 10000.3 PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Problémy konečné reprezentace ► uvažujme přesnost na 6 cifer ► součet 10000.0 + 3.14159 + 2.71828 ► přesně 10005.85987, zaokrouhleně 10005.9 naivně 10000.3 ► Kahanův algoritmus float input[] ; float sum = 0, err = 0; for (int i = 0; i= fabs(input[i])) err += (sum - t) + input[i]; el se err += (input[i] - t) + sum; sum = tmp; } err; ^^^^^^^^^^^^^^^|^^H^^^ PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 15/72 Problémy konečné reprezentace Robustní aritmetiky ► není to zadarmo ► Kahan 4 x, Neumaier 5x více operací ► třídění O (n log n) ► rekurzivní sčítání neeliminuje chybu zcela ► pozor na kompilátor: -funsafe-math-optimizations pro gcc PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Problémy konečné reprezentace Násobení a dělení ► samo o sobě nezanáší významnou chybu ► zachovává existující chybu PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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ší PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 ► chyba 0.082 %, tedy téměř o řád menší ► za cenu 6 aritmetických operací místo 2 ► bude 3x pomalejší < □ ► < igi ► < 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 ► chyba 0.082 %, tedy téměř o řád menší ► za cenu 6 aritmetických operací místo 2 ► bude 3x pomalejší ... uvidíme < □ ► < igi ► < Měření rychlosti výpočtu Naivní přístup ► POSIX volání gettimeofdayO gettimeofday(&start,NULL); c = a + b; c *= c; gettimeofday(&stop,NULL); ► současné CPU 2-4 GHz ► 1 aritmetická operace ~ 1-10 ns ► nemáme tak přesné hodiny PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Měření rychlosti výpočtu Zopakujeme v cyklu PA081: Programování numerických výpočtů A. Křenek ► opakování 109x navýší čas na měřitelnou ~ 1 s gettimeofday(&start,NULL); foř (i=0; i<1000000000; i++) c = a + b; c *= c; } gettimeofday(&stop,NULL); { O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Měření rychlosti výpočtu Zopakujeme v cyklu ► opakování 109x navýší čas na měřitelnou ~ 1 s gettimeofday(&start,NULL); foř (i=0; i<1000000000; i++) c = a + b; c *= c; } gettimeofday(&stop,NULL); { ► počítá podstatně rychleji ► trochu podezřelé PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Měření rychlosti výpočtu Zopakujeme v cyklu ► optimalizující kompilátor gcc -03 -funrol 1-loops ► 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říliš PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí = ► 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); foř (i=0; i<1000000000; i++) { c = a + b; c *= c; i f (nikdy) brzda(&nikdy,&a,&b,&c); } Měření rychlosti výpočtu Jak to dopadlo PA081: Programování numerických výpočtů A. Křenek ► Intel i7-2600 3.4 GHz vzorec čas Mflop/s Mcyklů/s (a + b)2 0.61 3278 1639 a2 + 2ab + b2 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 ► ► O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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ý. Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 < □ ► < igi ► 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 denormaliz ováné číslo O.OOOOOOmmm x ÍO"37 ► další výpočty nepřesné a navíc citelně pomalejší PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 25/72 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. Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování 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"); PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 27/72 Problémy konečné reprezentace Porovnávání nevinny vypočet float a=1.505, bl=1.315, b2=1.695, c=a+a, d=bl+b2; printf("%f %f %s\nM,c,d,c == d ? "jsou stejná" : "jsou ruzna"); ► má překvapivý výstup 3.010000 3.010000 jsou ruzna PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 27/72 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\nM,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 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 Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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\ x x y y PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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\ x x y y ► opatrný („miserly") přístup ► intervaly se i částečně překrývají => tvrdíme x = y ► jinak porovnáme x ^ y PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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\ x 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áváme x + ex > y - ey. PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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. Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Příklad Kvadratická rovnice ► rovnice tvaru ax2 + bx + c = 0 ► školní vzorec x+ = -b ± Vfo2 - 4ac 2a PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 31/72 Příklad Kvadratická rovnice ► rovnice tvaru ax2 + bx + c = 0 ► školní vzorec -b ± Vfo2 - 4ac x± =--- 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 Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Příklad Kvadratické rovnice ► konkrétně pro b = 2000, c = 1 ► „přesné" hodnoty ► VĎ= 1999.9989999997499998749999218749453... ► x+ = -.00050000012500006250003906252734377... ► x- = -1999.9994999998749999374999609374... ► pro float ► VĎ= 1999.999 ► x+ = -0.00048828125, chyba 2.5% ► X- = -1999.9995, chyba odpovídá přesnosti typu Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Příklad Kvadratické rovnice ► bezproblémové řešení pro x+ ► vlastnosti kořenů rovnice: ax2 + 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-3 vypočte se symetricky Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Příklad Kvadratické rovnice ► bezproblémové řešení pro x+ ► vlastnosti kořenů rovnice: ax2 + 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-3 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)deíinice ► „metoda počítá sice špatně, ale jenom trochu špatně" ► žádoucí a přitom realistická vlastnost všech numerických algoritmů PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Numerická stabilita (Pseudo)deíinice ► „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 ■X ■y PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Numerická stabilita (Pseudo)deíinice ► „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 ► algoritmus je stabilní, e* existuje-li pro každé x malé ex tak, že ey je malé ► „malé e" znamená zpravidla |e| < \5x ■y PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Numerická stabilita (Pseudo)deíinice ► „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 ► algoritmus je stabilní, e* existuje-li pro každé x malé ex tak, že ey je malé ► „malé e" znamená zpravidla |e| < \5x ■y PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí ► silnější definice zpětné stability - ey = 0 34/72 Numerická stabilita (Pseudo)deíinice ► „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 ► algoritmus je stabilní, e* existuje-li pro každé x malé ex tak, že ey je malé ► „malé e" znamená zpravidla |e| < \5x ■y PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí ► silnější definice zpětné stability - ey = 0 ► pro speciální druhy problémů jiné definice stability ► např. numerické integrování - „systém nevybuchne" 34/72 Numerická stabilita Příklady ► addss je numericky stabilní implementace sčítání y = xi + x2 ► zvolíme ô = 2~22 ► pro daná %\, y\ označme chybu jejich reprezentace 5i,2 ► ve float platí |5i,21 < 225 ► relativní chyba tohoto konkrétního sčítání je |ôa| < 2~25 y = adds(xi,x2) = (xi + x2)(l + ôa) = Xi(l + (5i)(l + 5a) + X2(l + 52)(1 + Sa) = (Xi + X2) + Xi((5i + Ôa + <5i(5a) + X2((52 + 5a + Ô2Ôa) ► i v nejhorším případě |5ij2 + ôa + <5i,25a| < 3 x 2~25 ► pro nejhorší případ x\ = x2 nastane \ý - y\ = \2xi(ôi + ôa + 5i5a)| < Numerická stabilita Příklady ► sčítání řady čísel 1230+1 + 1 + 1 + ...je nestabilní ► stabilní alternativa - setřídit a začít od nejmenšího PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 36/72 Numerická stabilita Příklady ► sčítání řady čísel 1230+1 + 1 + 1 + ...je nestabilní ► stabilní alternativa - setřídit a začít od nejmenšího ► úvodní příklad En = 1 - nEn-\ ► pro dostatečně velké N položíme EN = 0 ► počítáme v l-En %-l - - n ► v každém kroku je chyba redukována faktorem 1 In Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Zaokrouhlování ► žádoucí chování implementace aritmetických operací + ,-,*,/,^: op(x,y) = (x opy){\+5) pro |5| < ^(strojová přesnost) ► výsledek je stejně dobrý jako zaokrouhlený výsledek přesného výpočtu PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Zaokrouhlování ► žádoucí chování implementace aritmetických operací + ,-,*,/,^: op(x,y) = (x opy)(l+5) pro |5| < ^(strojová přesnost) ► výsledek je stejně dobrý jako zaokrouhlený výsledek přesného výpočtu ► stačí výpočet v dané přesnosti? PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Zaokrouhlování ► žádoucí chování implementace aritmetických operací + ,-,*,/,^: op(x,y) = (x opy)(l+5) pro |5| < ^(strojová přesnost) ► výsledek je stejně dobrý jako zaokrouhlený výsledek přesného výpočtu ► stačí výpočet v dané přesnosti? ► např. 0.100 x 21 - 0.111 x 2° (dvojkově, 3 cifry přesnosti) ► po zarovnání řádu: (0.100 - 0.011) x 21 = 0.001 x 21 ► správný výsledek byl 0.0001 x 21, chyba je 5 = 1 > 0.001 = u ► potřebujeme počítat s vyšší přesností ► (0.1000 - 0.0111) X 21 = 0.0001 X 21 ► tzv. guard digit, pro +, - stačí jedna PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Zaokrouhlování ► základní zaokrouhlení - k nejbližšímu reprezentovatelnému číslu 1.24 = 1.2 1.26 = 1.3 ► při nejednoznačnosti: ► vždy nahoru ► vždy dolů ► k sudému: 1.25 = 1.2, ale 1.35 = 1.4 ► asymetrické zaokrouhlení vnáší nežádoucí systematickou chybu ► např. průměr velkého souboru čísel vyšší/nižší Elementární funkce ► goniometrické, hyperbolické, exponenciální, logaritmické, odmocniny,... PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 39/72 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 PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Elementární funkce Co je špatně? ► celý aparát vede na iracionální čísla (rr, e, V2,...) ► 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 — y (a perioda) roste nade všechny meze ► sinx a cos x ► pro x — y 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 Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Dilema tvůrce tabulek ► počítáme např. ex na 3 místa, vychází (po n iteracích) 0.23450000... ► je to zaokrouhlením přesného výsledku 0.2344999999, správné zaokrouhlení 0.234 ► nebo jen nepočítám dostatečně přesně 0.23450001, správné zaokrouhlení 0.235 PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Dilema tvůrce tabulek ► počítáme např. ex na 3 místa, vychází (po n iteracích) 0.23450000... ► je to zaokrouhlením přesného výsledku 0.2344999999, správné zaokrouhlení 0.234 ► nebo jen nepočítám dostatečně přesně 0.23450001, správné zaokrouhlení 0.235 ► lze ukázat, že ex nemůže být strojově reprezentovatelné pro strojově reprezentovatelné x ► tedy k rozhodnutí stačí konečný počet kroků ► pro různé funkce různě, zpravidla cca. dvojnásobná přesnost PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí o 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ě ► čistě algebraické úpravy ► podobné (a + b)2 = a2 + 2ab + b2 ► podle kontextu zvolit vhodnou implementaci ► knihovní funkce ► vlastní (Taylorův rozvoj) uzpůsobené specifickým podmínkám PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Algebraické úpravy ► výraz Vl + x - Vl - x ► pro malá x odčítání velmi blízkých čísel ► transformace Vl + x - Vl - x = (VI + x - Vl - x)(Vl + x + Vl - x) 1 + x + VI - x (1 +x) - (1 -x) 1 + x + VI - x 2x 1 + x + VI - x PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí ► sčítání velmi blízkých čísel - dostatečně přesné Algebraické úpravy PA081: Programování numerických výpočtů A. Křenek ► výraz ln vx+T - ln ^fx ► pro velká x rozdíl blízkých čísel ► transformace ln Vx + 1 - ln Jx = = ln(x + 1)2- lnx? -(ln(x + 1) L x + 1 -ln-= 2 x - lnx) ^ln(l + -) 2 x ► možná ztráta přesnosti, ale i tak lepší O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 44/72 Mocninné řady Taylorův rozvoj PA081: Programování numerických výpočtů A. Křenek ► jedna ze základních vět matematické analýzy f(b) =f(a)+f'(a)^—^ + ,,, Ab-a)2 + f"(a) 2! 1! + + kde § je dostatečně malé (n) (b-a)n + O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Mocninné řady Taylorův rozvoj ► jedna ze základních vět matematické analýzy f(b)=f(a)+f'(a)^-^ + „ Ab-a)2 + f (a)---+ 2! in) (b ~ a)n + + kde § je dostatečně malé ► základ implementací většiny knihovních funkcí ► v extrémních případech potřebujeme „rozepsat" a provést algebraické úpravy PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 45/72 Mocninné řady Taylorův rozvoj základních funkcí ex = 1 + smx = x - cos x = 1 - ln(l + x) = x - X X2 2! X3 X5 3! + 5! X4 + 2! 4! X2 X3 + 2 3 + ... Vl + x = 1 + ^ - ^ 2 2-4 2-4-6 xi + 3x _L_ x 3x2 3 ■ 5x3 ITx ~ " 2" + 2^4 " 2 ■ 4 ■ 6 + PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 46/72 Mocninné řady Demonstrační příklad vyraz 1 - cosx x2 se pro x — 0 blíží i 2 ► 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 - cosx x2 = 0.47683709 ► 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 cos* = l-- + --- + ... ► pro výpočet na 8 cifer stačí do řádu x4 1 - COSX _l-l + -2"-24_l XZ x~2 " x~2 " 2 " 24 ► 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"? Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování 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 float a výsledky ~ 1 je to 10~7 ► při méně přesných vstupních datech méně striktní ► hledáme c aby pro x < c byl největší zanedbaný člen řady menší než požadovaná přesnost PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování 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 float a výsledky ~ 1 je to 10~7 ► při méně přesných vstupních datech méně striktní ► hledáme c aby pro x < c byl největší zanedbaný člen řady menší než požadovaná přesnost ► konkrétně v předchozím příkladu X < 1(T7 tedy x < 0.09211 6! ► 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.00000009997574124493 PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Mocninné řady ► výsledný kód pro O čem to bude float x,y; i f (x < 0.09211) y = 0.5 + x*x/24.0; el se y = (l-cos(x))/(x*x); Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 kde r = ||a - b ► silový vektor F = F a - b r ► snadný výpočet r2 = J^iat-bi) PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů ► r už vyžaduje odmocninu Shrnutí 52/72 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ě Vx 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 ^1+x ► iterační metody (Newtonova atd.) PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 53/72 Odmocnina Fast inverse square root ► hra Quake III na konci 90. let ► převratně realistická grafika ► vděčila velmi rychlému výpočtu Vx ► kód zřejmě pochází z grafických knihoven SGI i PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Odmocnina Fast inverse square root ► hra Quake III na konci 90. let ► převratně realistická grafika ► vděčila velmi rychlému výpočtu Vx ► kód zřejmě pochází z grafických knihoven SGI i float invSqrt(float x) { float xhalf = 0.5f*x,y; union { float x; i nt i ; } u; u.x = x; u.i = 0x5f3759df - (u.i » 1); y = u.x * (1.5f - xhalf * u.x * u.x); return y; } Odmocnina Fast inverse square root ► funguje na základě IEEE reprezentace čísla x = (1 + mx) x 2ex sign exponent (8 bits) _ significand (23 bits) I r 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) .23 0 ► pritom mx = Mx/2Zá a ex = Ex - 127 ► základ výpočtu 1 y = log2 y = - 2 loS2 * 1 log2(l + my) + £3, = --log2(l + mx) 1 2e* PA08i: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Odmocnina Fast inverse square root ► pro m e [0,1) si lze dovolit aproximaci log2(l +m) = m + a ► dosazením a úpravami dostaneme Ey x 223 +My = fl-|(£xx 223 +Mx) tj. červený řádek v kódu s empiricky stanovenou konstantou R ► více viz http://en.wikipedia.org/wiki/Fast i nverse_square_root Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Odmocnina Fast inverse square root PA081: Programování numerických výpočtů A. Křenek ► poslední krok je jedna iterace Newtonovy metody ► pro vstup x hledáme takové y aby y = ► tj. hledáme kořen funkce f(y) = yz-x ► Newtonovou metodou y-n+i = y n f (yn) f'(yn) = y yn n ^2 3 yn = yn + yn(l -xy2) což už odpovídá poslednímu výrazu v kódu O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Odmocnina Fast inverse square root ► aproximace je překvapivě přesná ► v Quake III je zakomentována další Newtonova iterace ► nebyla potřeba ► cca. 4x rychlejší než dělení ► dnes překonána instrukcí rsqrtss Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Odmocnina Fast inverse square root ► aproximace je překvapivě přesná ► v Quake III je zakomentována další Newtonova iterace ► nebyla potřeba ► cca. 4x rychlejší než dělení ► dnes překonána instrukcí rsqrtss ► přesné pochopení konkrétního účelu výpočtu ► včetně zhodnocení „má smysl tvrdě optimalizovat" ► volba adekvátní metody ► přispěla k velkému komerčnímu úspěchu ► pro jiné účely by byla zcela nevhodná Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Odmocnina Lze se jí i vyhnout ► Lenard-Jonesuv potenciál ► nevazebná interakce dvou atomů ► významná pro modelování biochemických dějů ► vzorec pro energii (potenciál) Elj = A B r 12 r 6 kde r je vzdálenost atomů 5 li 1 i i 4 - 3 - S 2 > - 1 - -e i l i ° r =21/6o 1-5 2 2.5 ) Odmocnina Lze se jí i vyhnout PA081: Programování numerických výpočtů A. Křenek ► výpočet velikosti síly F = dE dr 12 A 6B + r 13 r 7 ► liché exponenty => bude potřeba odmocnina O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Odmocnina Lze se jí i vyhnout PA081: Programování numerických výpočtů A. Křenek ► výpočet velikosti síly F = dE dr 12 A 6B yl3 y7 ► liché exponenty => bude potřeba odmocnina ► zajímá mě většinou silový vektor F = F-= (a - b) r r8 r14 ► tedy vystačím s r: O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 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 = r4 = 1.0/n; r2*r2; r8 = r4*r4; rl4 = r8*r4*r2; n = 12*A*rl4 - 6*B*r8; f or (i=0; i<3; F [i] *= n; Odmocnina Lze se jí i vyhnout ► implementace oat n,r2,r4,r8,rl4,F[4]; i nt i ; for (n=0,i=0; i<4; { F[i] = a[i] - b[i]; n += F[i]*F[i] ; } r2 r4 1.0/n; r2*r2; r8 = r4*r4; rl4 = r8*r4*r2; n = 12*A*rl4 - 6*B*r8; for (i=0; i<4; F [i] *= ► kvůli vektorizaci může mít smysl F [4] PA081: Programování numerických výpočtů A. Křenek O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí 61/72 Parametrizace rotace Aneb jak se vyhnout goniometrickým funkcím PA081: Programování numerických výpočtů A. Křenek ► ve 2D jeden úhel ► aplikace násobením maticí R = cos O čem to bude Čísla v plovoucí řádové čárce Numerická stabilita Zaokrouhlování Elementární funkce Algebraické úpravy Mocninné řady Odmocnina Algebra kvaternionů Shrnutí Parametrizace rotace Aneb jak se vyhnout goniometrickým funkcím ► ve 2D jeden úhel ► aplikace násobením maticí R = cos ► goniometrickým funkcím se můžeme vyhnout 2t 1 - t2 sin ch =--- cos ch =--- 1 l + t2 l + t2 ► numericky příjemné, nevyjádříme cj) = ±7 ► nelineární vztah