Transcript Cordic, dzielenie

Sprzętowa Implementacja
Algorytmów
Ernest Jamro, AGH Kraków
Cordic, Dzielenie, Pierwiastek...
Cordic (COordinate Rotation
DIgital Calculation)
Zobacz: http://my.execpc.com/~geezer/embed/cordic.htm
Omondi, A. R. „Computer arithmetic systems...”
Obrót o kąt 
X2 = X1 * cos() - Y1 * sin()
Y2 = X1 * sin() + Y1 * cos()
Po przekształceniu:
X2 = cos() * [ X1 - Y1 * tan() ]
Y2 = cos() * [ X1 * tan() + Y1 ]
2
CORDIC – c.d.
Wybieramy tg() tak aby = 1, ½, ¼, ...
Wtedy:
tan(21) = 1/1 21 = 45°
tan(32) = 1/2 32 = 26.5650°
tan(43) = 1/4 43 = 14.0362°
tan(54) = 1/8 54 = 7.12502°
tan(65) = 1/16 65 = 3.57633°
tan(76) = 1/32 76 = 1.78991°
cos(21) = 0.707107
cos(32) = 0.894427
cos(43) = 0.970142
cos(54) = 0.992278
cos(65) = 0.998053
cos(76) = 0.999512
W rezultacie:
Mnożenie przez tan()- zastępujemy przez dodawanie z
przesunięciem
Mnożenie przez cos() jest zastępowane przez pojedyncze
mnożenie skumulowane przez stałą wartość (zob. następna strona)3
CORDIC – Iteracje
Iteracja 1
X2 = cos( 21) * [ X1 - Y1 * tan( 21) ]
Y2 = cos( 21) * [ X1 * tan( 21) + Y1 ]
Iteracja 2
X3 = cos( 32) * { X2 - Y2 * tan( 32) }
Y3 = cos( 32) * { X2 * tan( 32) + Y2 }
itd.
Mnożenie przez cos() można wykonać na samym końcu a
nie w każdej iteracji.
4
CORDIC Przykład
Znaleźć sin(28.027°)
Stan początkowy:
 = 0° cos() = 1 X = 1
sin() = 0 Y = 0
Rotacja o 45°
X' = X - Y / 1 = 1 - 0 / 1 = 1
Y' = X / 1 + Y = 1 / 1 + 0 = 1
Rotacja z 45° do 18.435° ( 32 = -26.565°)
(Przesunięcie o kąt ujemny ale cos(- )= cos() więc nie ma to wpływu na wartość
skumulowanego iloczynu. Obracamy o kąt ujemny wtedy kiedy skumulowany kąt
obrotu jest większy od kąta docelowego.)
X' = X + Y / 2 = 1 + 1 / 2 = 1.5
Y' = -X / 2 + Y = -1 / 2 + 1 = 0.5
Rotacja z 18.435° do 32.471° ( 43 = 14.036°)
X' = X - Y / 4 = 1.5 - 0.5 / 4 = 1.375
Y' = X / 4 + Y = 1.5 / 4 + 0.5 = 0.875
5
CORDIC Przykład - c.d.
Rotacja z 32.471° do 25.346° ( 54 = -7.125°)
X' = X + Y / 8 = 1.375 + 0.875 / 8 = 1.484375
Y' = -X / 8 + Y = -1.375 / 8 + 0.875 = 0.703125
...
Rotacja 27.132° do 28.027° ( 87 = 0.895°)
X' = X - Y / 64 = 1.465300 - 0.750884 / 64 = 1.453567
Y' = X / 64 + Y = 1.465300 / 64 + 0.750884 = 0.773779
W rezultacie:
sin(28.027°) = 0.607253 * Y = 0.46988
cos(28.027°) = 0.607253 * X = 0.88268
Zamiast mnożyć przez skumulowany czynnik cos() można ustawić:
X= 0.607253 zamiast 1 w pierwszej iteracji.
Kąt  w każdej iteracji powinien być zapisany w pamięci LUT (w
przypadku implementacji sekwencyjnej) w przypadku implementacji
kaskadowej jest on wartością stałą.
6
Sprzętowa implementacja
y0
x0
z0
sign
>> 0
>> 0
b
b
a
const.
a
Sekwencyjna
sign
>> 1
>> 1
const.
y0
x0
sign
>> 2
FF
>> 2
z0
FF
FF
const.
sign
a
>> n
>> n
b
b
ROM
a
+/-
+/-
sign
>> 3
>> 3
const.
sign
>> 4
>> 4
const.
Kaskadowa
yn
xn
zn
7
Operacja dzielenia
(Algorytm: non-performing division)
mux
znak
01100100:1011
Dziesiętnie:
100:11= 9 reszta 1
+1 0 1 0 1
=
100= (01100100)2
00001
1
11= (01011)2
00011
-11= (10101)2
+1 0 1 0 1
9= (1001)2
<<1
11000
0
00110
Dzielna
+1 0 1 0 1
Dzielnik
11011
0
01100
+1 0 1 0 1
Odejmowanie
0 0 0 0 1=Reszta
znak
Wynik
rejestr przes.
1
8
Pozycjonowanie (normalizacja)
argumentów dzielenia
Aby przyspieszyć operacje dzielenia należy najpierw odpowiednio
przesunąć bitowo dzielną i dzielnik – tak aby na najbardziej
znaczących bitach występowały ‘1’
Przykład:
a) 1001: 0011 należy zastąpić przez 1001:1100 a wynik podzielić
przez 4 ( przesunąć w prawo o 2 bity) w ten sposób zyskujemy 2
takty zegara
b) 01100: 000101 należy zastąpić przez 11000 ( <<1) oraz 101000
(<<3) czyli wynik należy podzielić przez 4 (>>2)
Dokładniejszy opis normalizacja można znaleźć w wykładzie
liczby zmiennoprzecinkowe
9
Operacja dzielenia
(Przesunięcie po zerach)
01100100:1011
Dziesiętnie:
100:11= 9 reszta 1
+1 0 1 0 1
100= (01100100)2
11= (01011)2
-11= (10101)2
9= (1001)2
0001
Dwa zera wiec
przesuń dodatkowo
o 2 pozycje w lewo
i wpisz 2 zera do
wyniku
=
100
01100
10101
1
0 0 0 0 1 Reszta
10
Operacja dzielenia
Przesunięcie po zerach
<<1
Trzy najstarsze bity odejmowania (dzielnadzielnik):
<<2
010: wynik<<=1; wynik|=1; dzielna<<=1;
mux
001: wynik<<=2; wynik|= 2; dzielna<<=2;
MSBs
000: wynik<<=3; wynik|= 4; dzielna<<=3;
<<1
1xx: wynik<<=1; dzielna<<=1;
Dzielna
Dzielnik
Odejmowanie
MSBs
Przesunięcie o więcej niż jeden bit w
momencie kiedy wynik odejmowania
wskazuje, że Dzielnik > 2*Dzielna lub
Dzielnik> 4*Dzielna
Uwaga: argumenty muszą być
znormalizowane na samym początku
Wynik
rejestr przes.
11
Operacja dzielenia
Radix 4
Multiplekser:
If (w3>=0) then w= w3
Dzielnik
Else if (w2>=0) then w=w2
Dzielna
2
3
Else if (w1>=0) then w= w1
Else w= w0
<<2
w0
w1
w2
w3
mux
mux
w
12
Operacja dzielenia: Newton-Raphson
Znajdowanie miejsca zerowego funkcji:
f (Xi )
f '(Xi ) 
X i  X i 1
Nachylenie stycznej w
danym punkcie, f(Xi+1)=0
Przekształcając otrzymujemy
X0
X2
X1
f (Xi )
X i 1  X i 
f '( X i )
Szukamy miejsca zerowego funkcji:
1
f (X )  D 
X
Rozwiązanie (miejsce zerowe): X= 1/D
Pochodna f ’(X)= 1/X2
f (Xi )
D 1/ X i
X i 1  X i 
 Xi 
 X i (2  DX i )
2
f '(Xi )
1/ X i
13
Dzielenie: metoda Newtona – Raphsona
xi 1  xi (2  xi D)
Przykład: dzielenie przez 3, D=3
x0= 0.5
x1= 0.5(2 - 0.5·3)= 0.25
x2= 0.3125
x3= 0.33203125
x4= 0.333328247
x5= 0.333333333
Każda iteracja podwaja liczbę
poprawnych bitów wyniku
1
1
Ei 1   xi 1   (2  xi  D  xi2 ) 
D
D
 D(
1 2  xi
1

 xi2 )  D  (  xi ) 2 
D
D
D
 D Ei2
14
Newton-Raphson c.d.
:
xi 1  xi (2  xi b)
•Określić wartość wstępną x0 = 1/b np. za pomocą LUT i odpowiednich
przesunięć
•Iteracyjnie obliczyć xi+1 = xi(2-xib) aż do punktu zbieżności (w n
krokach). Proces iteracyjny jest tym szybciej zbieżny im wartość
początkowa x0 poszukiwanego rozwiązania jest bliższa wartości
rzeczywistej x.
Wymagana liczba iteracji, n, zależy od tego, z jaką precyzją obliczenia
mają zostać wykonane. Należy przy tym pamiętać, że każdy następny
krok iteracji podwaja ilość poprawnych bitów w wyniku. Można
założyć, że po czterech iteracjach jesteśmy w stanie obliczyć
odwrotność z dokładnością do 16 bitów. A po 5 iteracjach z
dokładnością 32 bitów.
15
Newton-Raphson c.d.
Dzielnik
xi 1  xi (2  xi b)
FF1
xib
2
LUT
FF – rejestry (ang.
Flip-Flops)
x0 = 1 / b
"sterowanie"
FF3
a
b
MUX1
a
b
MUX2
FF2
16
Pierwiastek kwadratowy
Newton-Raphson
f (Xi )
X i 1  X i 
f '( X i )
f (X )  X 2  P
f '(X )  2X
X i2  P 1
P
X i 1  X i 
 (Xi  )
2 Xi
2
Xi
Szybkość aproksymowania kwadratowa (liczba
ważnych bitów podwaja się w każdej iteracji)
Wada: - używanie operacji dzielenia
17
Pierwiastek kwadratowy bez
operacji dzielenia
Miejsce zerowe funkcji:
f '(X ) 
2
X3
1
f (X )  2  P
X
1
P
2
f (Xi)
Xi
X i 1  X i 
 Xi 
2
f '(Xi )
X i3
X
X i 1  i  (3  P  X i2 )
2
1
Wynikiem powyższej iteracji XK jest P dlatego aby otrzymać
bezpośrednio pierwiastek należy przeprowadzić operację dzielenia lub co
1
P
jest zalecane mnożenia:
P
 P X 
 P 18
XK
K
P
SQRT non-restoring digit-by-digit
19
Taylor-Maclaurin
f ( x)  f ( X  X )  f ( X ) 
f '(X )
f ''(X )
X 
X 2  ...
1!
2!
1
1 2
1 3 3
1 x  1 x 
x 
x  ..
2
24
246
2
3
x
x
e x  1  x    ...
2! 3!
x3 x5 x7
sin(x)  x     ...
3! 5! 7!
Ogólna postać: f(x)= a0+ a1x + a2x2 + a3x3 + ....
Postać ulepszona (nie trzeba obliczać xn) (metoda Horner’a):
f(x)= b0+ x·(b1 + x·(b2 + x·(b3 + ...) ) )
20
Różne metody wyznaczania
współczynników wielomianu
Wykres błędu aproksymacji dla różnych metod wyznaczania współczynników
wielomianu aproksymującego - dla funkcji ln(x), x=[1, 2) i wielomianu rzędu 3. Na
wykresie jest podany maksymalny błąd dla różnych metod
funkcja Matlab’a:
Minimax(Function, Interval, Degree)
Minimax – optymalizacja w celu osiągnięcia najmniejszego błędu w najgorszym przypadku
21
Wielomiany Czebyszewa
T0(x) =
1
T1(x) =
x
T2(x) =
2x2-1
T3(x) =
4x3-3x
T4(x) =
8x4-8x2+1
T5(x) =
16x5-20x3+5x
22
T6(x) =
32x6-48x4+18x2-1.
Redukcja dziedziny wejściowej
1) Redukcja zakresu danej wejściowej do mniejszego zakresu
2) Obliczanie funkcji dla zredukowanego przedziału
3) Rekonstrukcja funkcji dla całego przedziału wejściowego
Przykład:
Funkcja
przedział x
rekonstrukcja
sin(x)
[0, /2)
sin(x+k)= (-1)k sin(x), sin(-x)= sin(x)
2x
[0, 1)
2X= 2Xint2x
ln(x)
[1,2)
ln(X)= Xintln(2) +ln(x)
x
(0.25, 1]
(22ix)= 2i  x
1/x
[1, 2)
1/(2i x)= 2-i (1/x)
Xint- część całkowita X
23
Aproksymacja liniowa
Y=Yk+AK·(X-XK)
X
Zakłada się, że XK= MSB(X)
X-XK
XK
MSB
Y4
LUT
Y
LSB
LUT
A
YK
+
X1
X2
AK
*
X3
X4
Y3
YK=f(XK)
Y1
Y2
- interpolacja
AK=(Yk+1-Yk)/ 
funkcja skl ejana
(ang. spline)
= XK+1-XK
24
Aproksymacja kwadratowa
Y=Yk+BK·(X-XK) + AK(X-XK)2
X
X-XK LSB
XK
MSB
LUT
Y
LUT
B
BK
YK
X1
X2
Y2
Y1
+
LUT
A
*
AK
*
*
(X-XK)2
+
25
Aproksymacja kwadratowa 2
Inne rozwiązanie:
Y=YK + BK·(X-XK) + AK· (X-XK)·(XK+1-X)
= XK+1-XK
X=XK + x
Y=YK + BK·x + AK· x·(-x)
Wartości współczynników:
YK=f(XK)
BK=(Yk+1-Yk)/ 

4  f ( X K  )  2  f ( X K )  2  f ( X K 1 )
2
AK 
2
Założenie przy obliczaniu AK: f(XK+ /2)= Y(XK+ /2)
f(XK+ /2)= YK+BK ·/2+AK ·(/2)2
26
Aproksymacja kwadratowa 2
X
Y=YK + BK·x + AK· x·(-x)
Y= YK + x·[BK+ AK· (-x)]
x=X-XK
XK
LSB
MSB
LUT
Y
LUT
B
BK
YK
+
LUT
A
!x+1
AK
+
*
-x
*
27