Operace v FP a iterační algoritmy INP 2008 FIT VUT v Brně
1
Operace FP • Číslo X s pohyblivou řádovou čárkou X = MX.BEx zapíšeme jako dvojici (MX, EX), kde mantisa MX je ve dvojkovém doplňkovém kódu, nebo v přímém kódu se znaménkem na nM bitech, exponent EX je v kódu s posunutím 2nE-1 na nE bitech. Třetím definičním údajem je hodnota základu B. • Základní aritmetické operace na dvojici čísel X, Y s pohyblivou čárkou jsou: X + Y = (MX .2 Ex - Ey + MY).2 Ey , kde EX ≤ EY X - Y = (MX .2 Ex - Ey - MY).2 Ey , kde EX ≤ EY X * Y = (MX .MY).2 Ex + Ey X : Y = (MX : MY).2 Ex - Ey
2
Požadované operace v pohyblivé čárce • Z obecného zápisu je zřejmé, jaké operace jsou zapotřebí. • Pro sčítání a odčítání: – 1. Vypočte se v pevné čárce rozdíl EX - EY – 2. Posune se MX o EX - EY bitů (tj. doprava, pokud je EX ≤ EY) – 3. Vypočte se v pevné čárce MX .2 Ex - Ey +/- MY
• Dále je zapotřebí provést normalizaci a zaokrouhlení výsledku V = (MV, EV). Normalizovat výsledek znamená posouvat mantisu vlevo (vpravo) a podle toho zmenšovat (zvětšovat) exponent EV tak dlouho, až se do sledovaného bitu (v0 nebo v1) dostane 1. • Sčítání exponentů se provádí v binární sčítačce s korekcí, nebo ve speciální sčítačce pro posunutý kód s šířkou nE bitů, sčítání mantis se provádí v binární sčítačce se šířkou nM + 2 bitů s doplněným záchytným klopným obvodem S.
3
Princip obvodové realizace sčítačky/odčítačky FP EX
EY
Pokud EX ≤ EY, potom
MY
MX
1
2
X +/- Y = (MX .2 EX - EY +/- MY).2 EY
M s menším E M s větším E
7
1 3
SX SY
3
n = |EX - EY|
4 6
4
Obecný postup (nevíme, zda EX ≤ EY): 1 Rozdíl exponentů, nastavení sign a n 2 Prohození mantis podle sign
8
EX
3 Posuv mantisy o n bitů vpravo
EY
4 Sečtení/odečtení mantis podle Add/Sub
5
D
5 Zjištění počtu nul D v horních bitech mantisy 6
6 Normalizace mantisy, detekce spec. případů
7 7 Úprava exponentu podle D a sign
E’ - D SV
EV
MV
8 Logika pro výpočet znaménka 4
MX
Násobení čísel s pohyblivou řádovou čárkou
MY
EX
MV = MX * MY
EY
E V = EX + EY
Testování speciálních případů hodnot MV : MV = 0 MV nenormalizovaná MV normalizovaná
E V := 0
Posuv MV o 1 bit vlevo E V := E V - 1 Zaokrouhlení
X * Y = MX.MY 2 Ex+Ey Přeplnění?
výsledek V = (MV, EV)
ano
Posuv MV o 1 bit vpravo E V := E V + 1
ne Testování speciálních případů hodnotE V
EV > e
max
EV < e
min
:
_ e _ EV < e min < max
Konec
5
Obvodová realizace operací v FP • Pro operace s plovoucí čárkou prováděné v rámci uvedených algoritmů by bylo možno teoreticky použít sčítačku a násobičku ALU s pevnou čárkou. V praxi se to však takto nikdy neděje. Obvody aritmetiky s pohyblivou čárkou jsou konstruovány jako zcela nezávislá jednotka s vlastním řadičem.
6
Iterační algoritmy • Dělení čísel s pohyblivou čárkou a celou řadu dalších operací je možné v jednotce FPU provádět iteračními algoritmy. • Newtonův iterační algoritmus dělení – Tento algoritmus je založen na iteračním hledání průsečíku funkce f(x) s osou x pomocí tečny v bodě xi.
y
f(x)
xi
xi+1
x 7
Newtonův algoritmus dělení (1) Na základě odhadu xi hledáme přesnější odhad xi+1 v bodě průsečíku tečny funkce f s osou x. Rovnice přímky procházející bodem f(xi) je y - f(xi) = f '(xi)(x - xi) Pokládáme-li přímku za aproximaci funkce f(x), můžeme psát f(xi+1) - f(xi) = f '(xi)( xi+1 - xi) V průsečíku tečny s osou x je f(xi+1) = 0, takže odtud dostáváme iterační vzorec xi+1 = xi - f(xi)/f '(xi) Tento iterační vzorec nyní použijeme pro případ dělení. Máme-li dělit číslem b, zvolíme funkci f(x) = 1/x - b, a hledáme bod průchodu této funkce osou x, tedy bod x, kde f(x) = 0, takže pak 1/x = b, resp. x = 1/b. Operaci dělení číslem b pak nahradíme násobením číslem 1/b. Derivace f '(x) = -1/x2, odtud xi+1 = xi - (1/ xi - b)/(-1/ xi2) = = xi + xi - bxi2 = xi (2 - bxi)
8
Newtonův algoritmus dělení (2) Postup výpočtu a/b je následující: 1. Posune se b tak, aby padlo do intervalu (1, 2) a pomocí tabulky odhadů zvolíme první odhad x0 . 2. Provedeme krok iteračního výpočtu xi+1 = xi (2 - bxi). Krok 2 opakujeme tak dlouho, až se dosáhne požadované přesnosti na p bitů, kdy je relativní chyba (xi - 1/b)/(1/b) = 2-p. Počet správných (přesných) bitů se každým iteračním krokem zdvojnásobuje, tedy v následujícím kroku i+1 je relativní chyba ( xi+1 - 1/b)/ (1/b) = 2-2p 3. Výsledek n-té iterace xn vynásobíme číslem a, výsledek xn.a posuneme o odpovídající počet bitů podle kroku 1.
9
Newtonův algoritmus dělení - příklad Spočtěte binárně 1/b pro b = 20. (20)10 = (10100)2 xi+1 = xi (2 – bxi) • řádovou čárku posuneme tak, aby b ∈ <1, 2), tedy: 10100 → 1.0100 (o 4 bity doleva) • zvolíme např. x0 = 1, potom: x1 = x0(2 – bx0) = 1.(10 – 1,01) = 1.0,11 = 0,11 x2 = 0,11(10 – 1,01.0,11) = 0,11(10 – 0,1111) = 0,11.1,0001= 0,110011 x3 = 0,110011(10 – 1,01.0,110011) = 0,110011(10 – 0,11111111) = = 0,110011.1,00000001 = 0,11001100110011 atd. • řádovou čárku posuneme o 4 bity doleva: 0,000011001100110011 Ověření správnosti: 1/20 = 0,05 a (0,000011001100110011)2 = (0,051952362)10 10
Newtonův algoritmus dělení v HW xi+1 = xi (2 – bxi)
MUX - multiplexer REG_X– registr x REG_MINUS_B – registr -b MULT1 – násobička (-b*xi) ADDER– sčítačka (2,0+(-b)*xi) MULT2 - násobička (xi*(2,0+(-b)*xi)) FSM – konečný automat (kontrolér) SCNT – synchronní čítač
11
*
Goldschmidtův algoritmus dělení (1) Tento algoritmus se objevil poprvé asi u procesoru TI 8847, a prosadil se jako spolehlivý algoritmus dělení i dalších výpočtů s kvadratickou konvergencí. Byl použit v sálovém počítači IBM System/360, model 91, v jednotce FPU systému S/390, a v celé řadě navazujících serverů IBM, např. v serveru G4 v r. 1997 a dalších. Princip: zlomek a/b se rozšiřuje činitelem r zvoleným tak, aby se součin ve jmenovateli blížil k 1, tedy r .b → 1. Pak je r ≈ 1/b, a stačí vypočítat hodnotu součinu a.r . Rozšiřování zlomku: a r0 r1 . . . ri -1 b r0 r1 . . . ri -1
12
*
Goldschmidtův algoritmus dělení (2) Nyní položíme x0 = a, y0 = b. V iterační podobě počítáme v každém kroku výrazy xi +1 = ri xi, yi +1 = ri yi , kde xi a yi je i-tá aproximace hodnoty čitatele a jmenovatele po řadě. Z toho vyplývá, že podíl xi/yi = xi +1/yi +1 = konst. = a/b Když vybereme ri tak že yi → 1, pak xi → a/b, tedy xi se blíží k hledanému podílu. Pozn.: Stejný postup se dá použít k i výpočtu druhé odmocniny z a. Nechť x0 = a, y0 = a. V každém kroku vypočítáme xi+1 = ri2xi, yi+1 = ri yi. Pak xi+1/yi+12 = xi /yi2 = 1/a, takže když se ri vyberou tak, aby to vedlo na xi → 1, pak yi → √a. Tato metoda se používá k výpočtu druhé odmocniny v TI 8847.
13
*
Goldschmidtův algoritmus dělení (3) Jak zvolit ri? Víme, že x0 = a, y0 = b, xi +1 = ri xi a yi +1 = ri yi Předpokládejme, že se b blíží co nejvíce 1, ale je zatíženo chybou δ. Položme b = 1 - δ a IδI < 1. Zvolíme, že rozšiřující člen r0 se rovná r0 = 1 + δ ; protože y0 = b = 1 - δ , tak y1 = r0 y0 = (1 + δ) (1 - δ) = 1 – δ2 Opět položíme r1 = 1 + δ2, takže y2 = r1 y1 = 1 – δ4 a tak dále. Protože je IδI < 1, yi → 1. S tímto výběrem ri bude xi vypočítáno jako xi +1 = ri xi = (1 + δ2i) xi = [1 + (1 – b)2i ] xi, protože b = 1 - δ, odtud δ = 1 – b, což dosadíme do předchozí rovnice xi +1 = a [1 + (1 – b)] [1 + (1 – b)2] … [1 + (1 – b)2i ]
14
*
Goldschmidtův algoritmus dělení (upravený postup výpočtu pro rychlejší konvergenci)
Postup výpočtu: 1. a i b se posune tak, aby b leželo v intervalu <1, 2) 2. Vybere se odhad b* ≈ 1/b (z tabulky) 3. Položíme x0 = a b* , y0 = b b* (tímto trikem se zrychlí konvergence, protože bb* ≈ 1) 4. Iteruj, dokud xi není dostatečně blízko a/b. Cyklus: r ≈2–y
// aby platilo, že když je yi = 1 + δi, pak r ≈ 1 - δi
y = yr
// yi + 1 = yi r ≈ 1 - δi2
x =xr
// xi+1 = xi r
Konec cyklu.
0 1 2 3 4 5 6
x 9.00000000 4.5 5.625 5.9765625 5.99990845 6 6
y 1.50000000 0.75 0.9375 0.99609375 0.99998474 1 1
r x/y 0.50000000 6.00000000 1.25 1.0625 1.00390625 1.00001526 1 1
15
Další funkce Zde platí, že postupy vhodné pro programovou implementaci, jako MacLaurinův rozvoj nebo Čebyševovy polynomy atd. nemusí být pro obvodovou realizaci optimální, a proto byla odvozena řada modifikovaných nebo nových algoritmů. Druhá odmocnina Klasický postup hledání druhé odmocniny čísla A > 0 je založen na řešení rovnice x2 - A = 0. Aplikací Newtonovy iterační metody dostaneme xi+1 = 1/2 (A/xi +xi) Volíme x0 > 0 a po dostatečném počtu iterací dostaneme druhou odmocninu čísla A. Opět platí, že je-li i-tá iterace zatížena chybou d, tedy xi = √A(1 + d) pak v i+1 kroku dostaneme hodnotu xi+1 = √A(1 + d2/2(1+d)) Tedy každým krokem získáme dvojnásobný počet platných číslic.
16
Druhá odmocnina lépe Pro realizaci technickými prostředky je výhodnější hledat Newtonovou iterační metodou řešení rovnice 1/y2 - A = 0 Získáme tak iterační vzorec pro Newton-Raphsonův algoritmus yi+1 = yi (3 - A. yi2)/2 Volíme-li 0 < y0 < √ (3/A), pak po dostatečném počtu iterací dostaneme hodnotu 1/√A na požadovaný počet míst. Druhou odmocninu pak dostaneme výpočtem √A= A.(1/√A) (Druhou odmocninu lze počítat rovněž modifikovaným Goldschmidtovým algoritmem - TI 8847). Všeobecná zásada je najít co nejrychlejší algoritmy, založené pouze na operacích sčítání (odčítání), posuvů a případně i násobení.
17
CORDIC (1) Algoritmus CORDIC - Coordinate Rotational Digital Computer publikoval J. E. Volder v r. 1959. V poslední době se velmi rozšířil a byl značně modifikován pro další typy výpočtů. Jeho základní myšlenka umožňuje pomocí jednoho algoritmu počítat většinu matematických funkcí vyčíslováním funkce ve tvaru a ± b.2-i tedy pomocí součtů, rozdílů a posunů. Mezi nejznámější aplikace patří použití CORDICU v kapesní kalkulačce HP 35, v koprocesorech Intel počínaje I 8087, pro číslicovou Fourierovu transformaci, číslicovou filtraci signálů apod.
Metodu CORDIC můžeme označit za metodu skládání dílčích úhlů. 18
CORDIC (2) y B
yB r
y A
r
α
Φ
A
ϕ
V algoritmu Cordic je důležitý úhel Ø – např. jeho sinus a kosinus chceme vypočítat. Počáteční úhel ϕ je nastaven na nějakou konvenční hodnotu, např. 0. Přechod od ϕ k Ø se provádí v krocích – přičítání a odečítáním vhodných úhlů. Pokud vhodně zvolíme tyto úhly, stačí k výpočtu pouze sčítání, odčítání a posuny.
xB x x A Obecně vyjádříme přechod od bodu A k bodu B:
xA = r.cos ϕ yA = r.sin ϕ xB = r.cos (ϕ + α) = r.cos ϕ.cos α - r.sin ϕ.sin α yB = r.sin (ϕ + α) = r.sin ϕ.cos α + r.cos ϕ.sin α
(1) (2) 19
CORDIC (3) Dosadíme-li za r do (1) r = xA /cos ϕ, r = yA /sin ϕ pořadě, a do (2) v opačném pořadí, dostaneme pro cos α ≠ 0 xB = cos α (xA - yA .tg α) yB = cos α (yA + xA .tg α) Zvolíme-li xA = 1, yA = 0 (tj. ϕ = 0) a pro dané α určíme xB, yB (tedy otočíme vektor o α), dostaneme xB = cos α, yB = sin α. Požadovaný úhel natočení můžeme složit z n úhlů s kladným i záporným znaménkem (viz obrázek), tedy z orientovaných úhlů αi'. y y
Pozn. αi’ je buď +αi nebo -αi.
0
α2
α
3
α1 x =K 0
x 20
CORDIC (4) Výsledný úhel natočení α = α1'+ α2'+ ... αn' Postup výpočtu je pak následující: položíme x0' = xA, y0' = yA a určíme xi', yi' postupně pro i = 1, 2 ...n. Dostáváme pak iterační vzorce xi' = cos αi'.( xi-1' - yi-1' . tg αi') (3) yi' = cos αi'.( yi-1' + xi-1' . tg αi') (4) Ve dvojkové soustavě volíme takové úhly αi', že platí tg αi' = ± 21-i. Tím se v rovnicích (3), (4) nahradí násobení posuvem o i - 1 bitů. Další zjednodušení provedeme odložením násobení číslem cos αi' nakonec výpočtu, kdy výsledek vynásobíme číslem K = cos α1 . cos α2 ... cos αn kde αi = |αi' |, protože platí cos αi' = cos |αi' |. Pro i = 1, 2, ... n sestavíme tabulku hodnot αi, pro něž αi = arctg 21-i
21
Tabulka úhlů CORDIC i Tato tabulka se uloží do permanentní pamětí ROM.
0 1 2 3 4 5 6 7 8 9
α 45° 26,56 14,04 7,12 3,57 1,789 0,895 0,4476 0,2238 0,1119
tg α 1 0,5 0,25 0,125 0,0625 0,03125 0,015625 0,0078125 0,0039063 0,0019531
cos α 0,707 0,894 0,970 0,997 0,998 0,999 0,9998 0,9999 0,99999 0,999998
22
CORDIC (5) Pro x0 = 1, y0 = 0 tak dostaneme pro i = 1, 2, ...n konečný tvar iteračních vzorců xi = xi-1 - yi-1 . 21-i (5) yi = yi-1 + xi-1 . 21-i (6) Tyto vzorce platí, pokud bereme αi z tabulky s kladným znaménkem. Pokud má αi záporné znaménko, prohodíme v (5) a (6) operace plus a mínus. Pak cos α = K .xn , sin α = K .yn . Položíme-li x0 = K, y0 = 0, vyhneme se násobení a výsledek je xn = cos α , yn = sin α . Pro úhly větší než 90° využíváme symetrie a opakování grafu goniometrických funkcí. Příklad: viz cvičení
23
CORDIC – v jazyce C
y=sin(α)
http://en.wikipedia.org/wiki/CORDIC
int iterations; // počet iterací float arctanTable[iterations]; // tabulka úhlů v radiánech float K = 0.6073; // K float v_x=1,v_y=0; // složky vektoru v for(int i=0; i < iterations; i++) { arctanTable[i] = atan(pow(2,-i)); }
α x=cos(α)
// vytvoření tabulky Vstupní úhel je α (v
radiánech). Začíná se s vektorem v = (1,0).
// vlastní algoritmus Cordic for(i = 0; i < iterations; i++) { // pokud je α záporné, přičteme úhel z tabulky: if( α < 0) { v_x = v_x + v_y*pow(2,-i); v_y = v_y - v_x*pow(2,-i); α += arctanTable[i]; } // pokud je α kladné, odečteme úhel z tabulky else { v_x = v_x - v_y*pow(2,-i); v_y = v_y + v_x*pow(2,-i); α -= arctanTable[i]; } } v_x *= K; v_y *= K; 24
CORDIC – příklad výpočtu programu: Alfa = 0.523599 (rad) = 30.000001 (deg) Iteraci = 20, X = 1, Y = 0 Krok 0, Alfa = -0.261799, X = 1.000000, Y = 1.000000 Krok 1, Alfa = 0.201848, X = 1.500000, Y = 0.500000 Krok 2, Alfa = -0.043130, X = 1.375000, Y = 0.875000 Krok 3, Alfa = 0.081225, X = 1.484375, Y = 0.703125 Krok 4, Alfa = 0.018806, X = 1.440430, Y = 0.795898 Krok 5, Alfa = -0.012434, X = 1.415558, Y = 0.840912 Krok 6, Alfa = 0.003190, X = 1.428697, Y = 0.818794 Krok 7, Alfa = -0.004623, X = 1.422300, Y = 0.829955 Krok 8, Alfa = -0.000716, X = 1.425542, Y = 0.824400 Krok 9, Alfa = 0.001237, X = 1.427153, Y = 0.821615 Krok 10, Alfa = 0.000260, X = 1.426350, Y = 0.823009 Krok 11, Alfa = -0.000228, X = 1.425948, Y = 0.823705 Krok 12, Alfa = 0.000016, X = 1.426149, Y = 0.823357 Krok 13, Alfa = -0.000106, X = 1.426049, Y = 0.823531 Krok 14, Alfa = -0.000045, X = 1.426099, Y = 0.823444 Krok 15, Alfa = -0.000015, X = 1.426124, Y = 0.823401 Krok 16, Alfa = 0.000001, X = 1.426137, Y = 0.823379 Krok 17, Alfa = -0.000007, X = 1.426131, Y = 0.823390 Krok 18, Alfa = -0.000003, X = 1.426134, Y = 0.823385 Krok 19, Alfa = -0.000001, X = 1.426135, Y = 0.823382 X*K = 0.866092, Y*K = 0.500040 Výsledek cos(Alfa) = 0.866092, sin(Alfa) = 0.500040 Pozn.: 45 deg = 0,7854 rad a 26,56 deg = 0,4636 rad 25
CORDIC – v HW (základní implementace)
úhel v kroku i di je znaménko úhlu
26
CORDIC – další aplikace Metodu CORDIC lze použít též pro převod polárních souřadnic na pravoúhlé a naopak. Pro výpočet hodnot hyperbolických a hyperbolometrických funkcí byla vypracována modifikace algoritmu, kdy se koncový bod vektoru nepohybuje po kružnici, ale po hyperbole s rovnicí x2 - y2 = R2. Pro výpočet hodnot logaritmických funkcí se pak použije vztahu ln w = 2.arctgh (w - 1)/(w + 1) pro w > 1, případně ln w = 2.arccotgh (w + 1)/(w - 1) pro w < 1. Pro výpočet exponenciálních funkcí se použije vztahů ew = cosh w + sinh w aw = ew.ln a Cyklometrické funkce vyčíslujeme opačným postupem. Hodnotu arctg y/x, tedy úhel α pro x > 0 můžeme určit postupným natáčením vektoru (x, y) až do polohy, kdy y = 0. Při skládání úhlu α postupujeme opět podle vzorců (5), (6). 27
Porovnání časové náročnosti instrukcí (Intel, od 486, 487) Počty taktů při provedení instrukcí (1)
28
