Iterační algoritmy: Goldschmidtův pro dělení a CORDIC pro výpočet sin, cos Demonstrační cvičení 9 INP
Goldschmidtův iterační algoritmus pro dělení Pro výpočet a/b posuneme a, b na a', b' tak, že b' ∈ <1, 2). Položíme x0=a', y0=b' a rozšiřujeme zlomek xi/yi, zlomkem ri/ri po n kroků dokud není yn→1. Vyjádříme-li jmenovatel ve tvaru y0=1-δ, kde |δ|<1, pak v prvním iteračním kroku provedeme rozšíření: x1/y1 = (x0/y0)*(r0/r0) = [x0 / (1-δ)]*[(1+δ) / (1+δ)] = x0(1+δ)/(1-δ2) V dalších iteracích obecně: xi+1/yi+1 = (xi/yi) * (ri/ri) = (xi/yi) * [(1+δ2i)/(1+δ2i)] Můžeme též psát: yi = 1-δi a yi+1 = 1-δi2 Protože yi=1-δi a ri=1+δi, lze odvodit: ri = 2-yi
Goldschmidt: příklad dělení 16 / 5 = 3.2 A = 16(10) = 10000(2) b = 5(10 )= 101(2) a' = 100(2) = 4(10) b' = 1.01(2) = 1.25(10) = 1+δ = 1+0.25; δ=0.01(2) = 0.25(10) x0 = a ' y0 = b' r0 = 2-y0 = 2-1.25 = 0.75 x1 = x0r0 = 4 * 0.75 = 3 y1 = y0r0 = 1.25 * 0.75 = 0.9375 r1 = 2-y1 = 2-0.9375 = 1.0625 x2 = x1r1 = 3 * 1.0625 = 3.1875 y2 = y1r1 = 0.9375 * 1.0625 = 0.99609375 r2 = 2-y2 = 2-0.99609375 = 1.00390625 x3 = x2r2 = 3.1875 * 1.00390625 = 3.19995117188 y3 = y2r2 = 0.99609375 * 1.00390625 = 0.99998474121 r3 = 2-y3 = 2-0.99998474121 = 1.00001525879 x4 = x3r3 = 3.19995117188 * 1.00001525879 = 3.19999999926 y4 = y3r3 = 0.99998474121 * 1.00001525879 = 0.99999999976 r4 = 2-y4 = 2-0.99999999976 = 1.00000000023
CORDIC (Coordinate Rotational Digital Computer – J. E. Volder, 1959) Základní myšlenka: Umět vypočítat většinu matematických funkcí vyčíslováním funkce ve tvaru a±b.2-i tj. pomocí součtů, rozdílů, posunů, případně vyhledávání v look-up tabulce.
CORDIC v úloze výpočtu sin(Φ) a cos(Φ) S využitím vztahů pro výpočet sin(ϕ ϕ + α) a cos(ϕ ϕ + α) lze odvodit vzorce pro výpočet [xB, yB] ze známých souřadnic [xA, yA]: xB = xA * cos(α) - yA * sin(α) yB = xA * sin(α) + yA * cos(α)
[xB, yB]
V algoritmu CORDIC je důležitý úhel Φ – 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. Pokud vhodně zvolíme tyto kroky (tj. velikosti úhlu α v každém kroku), stačí k výpočtu pouze sčítání a rotace. Předchozí vztahy mohou být přepsány následovně: [xA, yA]
xB = cos(α) * [ xA - yA * tan(α) ] yB = cos(α) * [ yA + xA * tan(α) ]
CORDIC (2) xB = cos(α) * [ xA - yA * tan(α) ] yB = cos(α) * [ yA + xA * tan(α) ] Hodnoty pro α se volí tak, že tan(α) nabývá hodnot 1/(mocina čísla 2): i tg α cos α α 0
45°
1
0,707
1
26,56
0.5
0,894
2
14,04
0.25
0,970
3
7,12
0.125
0,997
4
3,57
0.0625
0,998
5
1,789
0,03125
0,999
6
0,895
0,015625
0,9998
7
0,4476
0,0078125
0,9999
To nám dovoluje nahradit násobení tan(α) rychlou operací posun vpravo.
Je však ještě třeba vypořádat se s cos(α α). Ukažme si několik po sobě jdoucích iterací:
CORDIC (3)
První iterace (z [X0, Y0] do [X1, Y1]) znamená otočení o úhel α0 X1 = cos(α0) * ( X0 – Y0 * tan(α0) ) xi+1 = cos(α) * [ xi – (yi >> i) ] Y1 = cos(α0) * ( Y0 + X0 * tan(α0) ) yi+1 = cos(α) * [ yi + (xi >> i) ] Druhá iterace (z [X1, Y1] do [X2, Y2]) znamená otočení o úhel α1 X2 = cos(α1) * ( X1 – Y1 * tan(α1) ) Y2 = cos(α1) * ( Y1 + X1 * tan(α1) )
Dosazením získáme: X2 = cos(α1) * { cos(α0) * [ X0 – Y0 * tan(α0) ] - cos(α0) * [ Y0 + X0 * tan(α0)] * tan(α1) } = cos(α1) * cos(α0) * { [ X0 – Y0 * tan(α0) ] - [ Y0 + X0 * tan(α0)] * tan(α1) } Je patrné, že kosiny se před závorkami násobí, tj. cos(α0)*cos(α1)*cos(α2)*...*cos(αn) Vyjádřením hodnot α jako inverzních hodnot tangenty dostaneme řadu, ze které potom hodnotu:
Tato hodnota se nazývá agregační konstanta. Po vynásobení agregační konstantou získáme požadovaný sinus a kosinus:
Př. S využitím CORDIC vypočtěte sin(28.027°) a cos(28.027°). Položíme: ϕ = β0 = 0° cos(ϕ) = 1 X0 = 1 sin(ϕ) = 0 Y0 = 0
xi+1 = xi – (yi >> i) yi+1 = yi + (xi >> i)
Přičti úhel α0 = 45° → β1 = β0 + α0 = 0 + 45 = 45° X1 = X0 – Y0 / 1 =1-0/1 =1 Y1 = Y0 + X0 / 1 =0+1/1 =1 Odečti úhel α1 = 26.565° → β2 = β1 - α1 = 45 - 26.565 = 18.435° Protože se jedná o záporný úhel a protože je tangens funkce lichá, změníme znaménko u posouvaných čísel: X2 = X1 + Y1 / 2 =1+1/2 = 1.5 Y2 = Y1 – X1 / 2 =1-1/2 = 0.5 Přičti úhel α2 = 14.036° → β3 = β2 + α2 = 18.43 + 14.036 = 32.471° X3 = X2 – Y2 / 4 = 1.5 - 0.5 / 4 = 1.375 Y3 = Y2 + X2 / 4 = 0.5 + 1.5 / 4 = 0.875 Odečti úhel α3 = 7.125° → β4 = β3 - α3 = 32.471 - 7.125 = 25.346° X4 = X3 + Y3 / 8 = 1.375 + 0.875 / 8 = 1.484375 Y4 = Y3 – X3 / 8 = 0.875 - 1.375 / 8 = 0.703125
Přičti úhel α4 = 3.576° → β5 = β4 + α4 = 25.346 + 3.576 = 28.922° X5 = X4 – Y4 / 16 = 1.484375 - 0.703125 / 16 = 1.440429 Y5 = Y4 + X4 / 16 = 0.703125 + 1.484375 / 16 = 0.795898 Odečti úhel α5 = 1.790° → β6 = β5 - α5 = 28.922 - 1.790 = 27.132° X6 = X5 + Y5 / 32 = 1.440429 + 0.795898 / 32 = 1.465300 Y6 = Y5 – X5 / 32 = 0.795898 - 1.440429 / 32 = 0.750884 Přičti úhel z α6 = 0.895° → β7 = β6 + α6 = 27.132 + 0.895 = 28.027° X7 = X6 – Y6 / 64 = 1.465300 - 0.750884 / 64 = 1.453567 Y7 = Y6 + X6 / 64 = 0.750884 + 1.465300 / 64 = 0.773779 Předpokládejme, že nám přesnost výsledku již postačuje a výpočet ukončíme. Následuje násobení agregační konstantou a získání požadovaných hodnot: sin(28.027°) = 0.607253 * Y = 0.46988 cos(28.027°) = 0.607253 * X = 0.88268 Abychom se vyhnuli tomuto násobení, můžeme počáteční hodnotu X nastavit na 0.607253 místo 1. Znamená to vlastně, že budeme otáčet vektor (mezi osami X a Y) o délce 0.607253 místo o délce 1. Při každé iteraci je třeba rozhodnout, zda přičíst nebo odečíst následující hodnotu α. To se provede porovnáním βn s cílovou hodnotou.
const int iterations = 22; float alpha; // angle in Radians float arctanTable[iterations]; // in Radians float K = 0.607253; float v_x = 1.0, v_y = 0.0; // Vector v; x and y components int i;
CORDIC v C
for (i=0; i < iterations; i++) arctanTable[i] = atan(pow(2, -i)); float vnew_x; // To store the new value of x for (i = 0; i < iterations; i++) { // If alpha is negative => do a counter-clockwise rotation if ( alpha < 0) { vnew_x = v_x + v_y * pow(2, -i); v_y = v_y - v_x * pow(2, -i); alpha = alpha + arctanTable[i]; } // If alpha is positive, we need to do a clockwise rotation else { vnew_x = v_x - v_y * pow(2, -i); v_y = v_y + v_x * pow(2, -i); alpha = alpha - arctanTable[i]; } v_x = vnew_x; } v_x = v_x * K; v_y = v_y * K;
CORDIC – příklad výpočtu programu alpha = 45.000000 deg. = 0.785398 rad. Iter. 0, alpha = 0.000000, x = 1.000000, y = 1.000000 Iter. 1, alpha = -0.463648, x = 0.500000, y = 1.500000 Iter. 2, alpha = -0.218669, x = 0.875000, y = 1.375000 Iter. 3, alpha = -0.094314, x = 1.046875, y = 1.265625 Iter. 4, alpha = -0.031895, x = 1.125977, y = 1.200195 Iter. 5, alpha = -0.000655, x = 1.163483, y = 1.165009 Iter. 6, alpha = 0.014968, x = 1.181686, y = 1.146829 Iter. 7, alpha = 0.007156, x = 1.172726, y = 1.156061 Iter. 8, alpha = 0.003250, x = 1.168210, y = 1.160642 Iter. 9, alpha = 0.001297, x = 1.165944, y = 1.162924 Iter. 10, alpha = 0.000320, x = 1.164808, y = 1.164062 Iter. 11, alpha = -0.000168, x = 1.164239, y = 1.164631 Iter. 12, alpha = 0.000076, x = 1.164524, y = 1.164347 Iter. 13, alpha = -0.000046, x = 1.164382, y = 1.164489 Iter. 14, alpha = 0.000015, x = 1.164453, y = 1.164418 Iter. 15, alpha = -0.000015, x = 1.164417, y = 1.164453 Iter. 16, alpha = -0.000000, x = 1.164435, y = 1.164436 Iter. 17, alpha = 0.000007, x = 1.164444, y = 1.164427 Iter. 18, alpha = 0.000004, x = 1.164439, y = 1.164431 Iter. 19, alpha = 0.000002, x = 1.164437, y = 1.164433 Iter. 20, alpha = 0.000001, x = 1.164436, y = 1.164434 Iter. 21, alpha = 0.000000, x = 1.164436, y = 1.164435 sin(45.000000) = 0.707107, cos(45.000000) = 0.707107
CORDIC v HW
xi+1=xi + yi.di.2-i yi+1=yi - xi.di.2-i
zi+1=zi-di.arctanTab[i] di=-1 if zi<0, +1 otherwise
