Últimos Cambios
Blog personal: El hilo del laberinto
|
|
|
Últimos Cambios |
|
Blog personal: El hilo del laberinto
|
|
|
Última Actualización: 26 de Febrero de 2003 - Miércoles
Tradicionalmente el cálculo de funciones transcendentales y trigonométricas se ha realizado mediante expansiones polinómicas, como polinomos de Taylor, Min-Max o Chevichev. Dejando a un lado sus características de convergencia y error, la evaluación de polinomios, aunque se simplifique mediante la regla de Horner, requiere el uso de multiplicaciones. En la actualidad, las CPUs modernas disponen de instrucciones de multiplicación muy rápidas, pero no es el caso en entornos empotrados o en desarrollos ASIC.
El algoritmo CORDIC (Coordinate Rotation Digital Computer), propuesto en 1959, permite el cálculo de este tipo de funciones de una forma muy simple, rápida y precisa, sin requerir hardware caro como una unidad de multiplicación.
La formulación normal de CORDIC, de forma genérica, es:
x[n+1]=x[n]-m*s[n]*y[n]*2-n y[n+1]=y[n]+s[n]*x[n]*2-n z[n+1]=z[n]-s[n]*e[n]
Donde m es 0, +1 o -1, s[n] es +1 o -1, y e[n] son constantes precalculadas.
Supongamos un vector bidimensional, con origen (0,0) y extremo (x,y). Si giramos el vector un ángulo A, respecto al origen (0,0), su nuevo extremo será (x',y')=(x*cos(A)-y*sin(A),x*sin(A)+y*cos(A)). Este nuevo vector, por supuesto, puede rotarse un ángulo adicional aplicando la fórmula las veces que sea necesario.
Si partimos de un vector (1,0) y lo giramos un ángulo A para obtener un vector (x,y), sabemos por trigonometría básica que cos(A)=x y sin(A)=y. Es decir no solo obtenemos senos y cosenos, sino que los podemos obtener simultaneamente.
Pero la ecuación propuesta no solo requiere multiplicaciones, sino que requiere el cálculo de esos senos y cosenos que son lo que queremos obtener. ¿Dónde está el truco?.
Partiendo de la ecuación, sabemos que para rotar un vector un ángulo A, tenemos, de forma iterativa: (x[n+1],y[n+1])=(x[n]*cos(A[n])-y[n]*sin(A[n]),x[n]*sin(A[n])+y[n]*cos(A[n])), donde (x[n],y[n]) es la serie de vectores que vamos obteniendo, y A[n] es la serie de ángulos que vamos aplicando a cada vector.
Si sacamos el coseno como factor común, tenemos (x[n+1],y[n+1])=cos(A[n])*(x[n]-y[n]*tan(A[n]),x[n]*tan(A[n])+y[n])).
El ángulo final A será la suma de ángulos A[n], aplicados etapa a etapa. Los valores concretos de cada A[n] son irrelevantes, y los podemos elegir como queramos, siempre que cubran el rango y la precisión que necesitamos. Tanto es así que podemos elegir ángulos A[n] que cumplan una propiedad que es la "idea feliz" que hace posible CORDIC.
Supongamos que elegimos los A[n] de forma que tan(A[n])=2-n. Con esto, el cálculo de la tangente desaparece, y la multiplicación con ella se reduce a un desplazamiento de bits, que es una operación tremendamente eficiente.
Completando las iteraciones y sacando factor común, obtenemos:
(cos(A),sin(A)) = (x,y) = Producto(cos(A[n]))*Iteración(x[n]-y[n]>>n),x[n]>>n+y[n]). El vector inicial será (1,0), que se corresponde a un ángulo de cero grados.
A la hora de calcular el ángulo A a partir de los A[n] existen dos opciones básicas. Supongamos que S[n] es el ángulo que llevamos rotado hasta el momento n. La primera opción es
S[n+1]=S[n]+A[n] si S[n+1] La segunda opción consiste en que S[n+1]=S[n]+M[n]*A[n],
donde M[n] puede ser el valor +1 o -1. La
elección de uno u otro puede hacerse comparando el error respecto al
ángulo final A. Otra opción simple es usar +1 si
estamos por debajo del ángulo final, y -1 si estamos por encima.
Dado que cos(A)=cos(-A), con esta elección
Producto(cos(A[n])) es una constante independiente del ángulo
final a obtener. Sólo depende de la elección de los A[n]. Valores que,
por otra parte, están ya elegidos.
Veamos un ejemplo en lenguaje Python (escrito para que se entienda, no para
que sea eficiente):
En el constructor se crean las tablas que vamos a necesitar. En una
implementación hardware, esas tables estaría ya integradas en el diseño
electrónico del sistema.
El método sincos() nos devuelve, simultaneamente, el seno
y el coseno del ángulo solicitado.
La precisión depende de:
Cuando se implementa por software en una CPU moderna, el mayor
problema es el tercer punto, ya que la mayoría de los algoritmos
de cálculo de seno y coseno (en particular IEEE-754) garantizan un
error bajo (con un gran coste computacional, eso sí). En ese
sentido una implementación CORDIC no puede competir en términos
de precisión, aunque el resultado es sorprendentemente bueno.
En mi sistema, la clase Python anterior tiene una precisión de cinco
dígitos decimales, en el peor de los casos. Esta precisión se puede incrementar
siendo consciente de las pérdidas de precisión de los números en coma
flotante cuando se restan números, por ejemplo, de órdenes de magnitud
muy diferentes. Pero éste es otro tema.
Si la implementación es hardware, el algoritmo CORDIC arrasa.
Este algoritmo se puede usar para rotar vectores arbitrarios, convertir
coordenadas entre notación polar y cartesiana (en ambos sentidos) y obtener
la funciones trigonométricas inversas (dado un coseno, por ejemplo, obtener
su ángulo). Modificándolo de forma obvia (y haciendo m=-1,
se puede usar la misma técnica para trabajar con trigonometría hiperbólica.
Si en la fórmula canónica de CORDIC hacemos m=1,
e[n]=tan-1(2-n),
s[n]=-sign(y[n]),
z[0]=0, obtenemos que:
z[n+1]=tan-1(y[0]/x[0])
y
x[n+1]=sqrt(x[0]2+y[0]2)/K, donde
K es el producto de los cosenos de e[n].
La rutina en Python sería:
Si en la fórmula canónica de CORDIC hacemos m=-1,
e[n]=tanh-1(2-n),
s[n]=sign(z[n]),
y[0]=0, x[0]=1 y z[0]=A,
obtenemos el seno y coseno hiperbólicos del ángulo A.
Para asegurar la convergencia, es necesario repetir algunos valores
de e[n]. En concreto hace falta repetir los índices de la forma
4, 13, 39, k, 3*k+1,...
La rutina en Python sería:
Si en la fórmula canónica de CORDIC hacemos m=-1,
e[n]=tanh-1(2-n),
s[n]=-sign(y[n]),
z[0]=0, obtenemos que:
z[n+1]=tanh-1(y[0]/x[0])
y
x[n+1]=sqrt(x[0]2-y[0]2)/K, donde
K es el producto de los cosenos hiperbólicos de e[n].
Para asegurar la convergencia, es necesario repetir algunos valores
de e[n]. En concreto hace falta repetir los índices de la forma
4, 13, 39, k, 3*k+1,...
La rutina en Python sería:
Si en la fórmula canónica de CORDIC hacemos m=-1,
e[n]=tanh-1(2-n),
s[n]=-sign(y[n]),
x[0]=w+1, y[0]=w-1 y z[0]=0,
obtenemos el logaritmo neperiano del valor w.
Para asegurar la convergencia, es necesario repetir algunos valores
de e[n]. En concreto hace falta repetir los índices de la forma
4, 13, 39, k, 3*k+1,...
La rutina en Python sería:
La velocidad de convergencia (error y número de pasos) depende
del valor concreto a calcular. Llegado el caso, resulta muy conveniente
realizar preescalados para aplicar el algoritmo en los rangos más
favorables.
Tradicionalmente las raíces cuadradras se calculan como
sqrt(x)=exp(log(x)/2). Por supuesto, esta opción
es tremendamente lenta.
Se puede implementar un sistema mucho más
eficiente, aunque requiere multiplicaciones.
Esta rutina obtiene la raíz cuadrada de un número inferior a
65536. El grado de precisión es prácticamente ilimitado, y depende
de la precisión de los cálculos intermedios y del número de
iteraciones que elijamos. Hay que tener mucho cuidado con el tema
de la precisión intermedia de los resultados en coma flotante.
Si tenemos valores muy grandes o muy pequeños, se usa un escalado,
y solucionado. Por ejemplo, si queremos sacar la raíz cuadrada
de 1224744871, lo dividimos primero por 65536, calculamos la raíz
cuadrada de ese valor y multiplicamos el resultado por 256.
El algoritmo CORDIC no requiere multiplicaciones.
Si hacemos m=-1,
e[n]=tanh-1(2-n),
s[n]=-sign(y[n]),
x[0]=w+0.25, y[0]=w-0.25 y z[0]=0,
obtenemos la raíz cuadrada de w.
Para asegurar la convergencia, es necesario repetir algunos valores
de e[n]. En concreto hace falta repetir los índices de la forma
4, 13, 39, k, 3*k+1,...
La rutina en Python sería:
La velocidad de convergencia (error y número de pasos) depende
del valor concreto a calcular. Llegado el caso, resulta muy conveniente
realizar preescalados para aplicar el algoritmo en los rangos más
favorables.
La rutina en Python sería:
Aunque la mayoría de las CPUs modernas incluyen estas funciones
como primitivas relativamente rápidas (la división suele seguir
siendo bastante lenta comparada con el resto de operaciones),
en entornos integrados o en CPUs pequeñas (PIC, microcontroladores)
se requiere un sistema simple y eficiente.
Podemos tener y[n+1]=x[0]*z[0] haciendo, en
la fórmula canónica de CORDIC, m=0,
e[n]=2-n, s[n]=sign(z[n]),
y[0]=0.
Para obtener z[n+1]=y[0]/x[0], hacemos
m=0, e[n]=2-n, s[n]=-sign(y[n]),
z[0]=0.
Los algoritmos resultantes son muy similares a los normales
de lápiz y papel.
Por ejemplo, a la hora de calcular el seno y el coseno, podemos obviar
las dos multiplicaciones finales si en vez de inicializar el sistema
con el vector (1,0), se inicializa con el vector
(self.c,0).
Más información sobre los OpenBadges
Donación BitCoin: 19niBN42ac2pqDQFx6GJZxry2JQSFvwAfS
class cordic :
def __init__(self) :
import math
angulos=[]
c=1.0
p=1.0
for i in xrange(16) :
a=math.atan(p)
c*=math.cos(a)
p/=2.0
angulos.append(a)
self.c=c
self.angulos=angulos
def sincos(self,angulo) :
x,y=(1.0,0.0)
p=1.0
a=0
for i in self.angulos :
if a<angulo :
(x,y)=(x-y/p,y+x/p)
a+=i
else :
(x,y)=(x+y/p,y-x/p)
a-=i
p*=2
return (y*self.c,x*self.c)
Arcotangente y magnitud
class cordic :
def __init__(self) :
import math
angulos=[]
c=1.0
p=1.0
for i in xrange(16) :
a=math.atan(p)
c*=math.cos(a)
p/=2.0
angulos.append(a)
self.c=c
self.angulos=angulos
def tan1mag(self,x,y) :
z=0.0
p=1.0
for i in self.angulos :
if y>=0 :
x,y,z=x+y/p,y-x/p,z+i
else :
x,y,z=x-y/p,y+x/p,z-i
p*=2
return z,x*self.c
Seno y coseno hiperbólicos
class cordic :
def __init__(self) :
import math,cmath
angulos=[]
c=1.0
p=0.5
for i in xrange(16) :
a=cmath.atanh(p).real
c*=math.cosh(a)
p/=2.0
angulos.append(a)
self.c=c
self.angulos=angulos
def sincosh(self,a) :
import math
x=1.0
y=0.0
z=float(a)
p=2
j=1
d=4
c=self.c
for i in self.angulos :
if z<0 :
x,y,z=x-y/p,y-x/p,z+i
else :
x,y,z=x+y/p,y+x/p,z-i
if d==j :
c*=math.cosh(i)
d=3*d+1
if z<0 :
x,y,z=x-y/p,y-x/p,z+i
else :
x,y,z=x+y/p,y+x/p,z-i
j+=1
p*=2
return y*c,x*c
Arcotangente hiperbólica y raíz cuadrada de la resta de cuadrados
class cordic :
def __init__(self) :
import math,cmath
angulos=[]
c=1.0
p=0.5
for i in xrange(16) :
a=cmath.atanh(p).real
c*=math.cosh(a)
p/=2.0
angulos.append(a)
self.c=c
self.angulos=angulos
def atanhmagneg(self,x,y) :
import math
x=float(x)
y=float(y)
z=0.0
p=2
j=1
d=4
c=self.c
for i in self.angulos :
if y>=0 :
x,y,z=x-y/p,y-x/p,z+i
else :
x,y,z=x+y/p,y+x/p,z-i
if d==j :
c*=math.cosh(i)
d=3*d+1
if y>=0 :
x,y,z=x-y/p,y-x/p,z+i
else :
x,y,z=x+y/p,y+x/p,z-i
j+=1
p*=2
return z,x*c
Logaritmos neperianos
class cordic :
def __init__(self) :
import math,cmath
angulos=[]
c=1.0
p=0.5
for i in xrange(16) :
a=cmath.atanh(p).real
c*=math.cosh(a)
p/=2.0
angulos.append(a)
self.c=c
self.angulos=angulos
def log(self,w) :
import math
x=float(w)+1
y=float(w)-1
z=0.0
p=2
j=1
d=4
c=self.c
for i in self.angulos :
if y>=0 :
x,y,z=x-y/p,y-x/p,z+i
else :
x,y,z=x+y/p,y+x/p,z-i
if d==j :
c*=math.cosh(i)
d=3*d+1
if y>=0 :
x,y,z=x-y/p,y-x/p,z+i
else :
x,y,z=x+y/p,y+x/p,z-i
j+=1
p*=2
return 2*z
Raíz cuadrada
def sqrt(valor) :
base=256.0
y=yy=0.0
for i in xrange(16) :
yy=y+base
if yy*yy <= valor :
y=yy
base/=2.0
return y
class cordic :
def __init__(self) :
import math,cmath
angulos=[]
c=1.0
p=0.5
for i in xrange(16) :
a=cmath.atanh(p).real
c*=math.cosh(a)
p/=2.0
angulos.append(a)
self.c=c
self.angulos=angulos
def sqrt(self,w) :
import math
x=float(w)+0.25
y=float(w)-0.25
z=0.0
p=2
j=1
d=4
c=self.c
for i in self.angulos :
if y>=0 :
x,y,z=x-y/p,y-x/p,z+i
else :
x,y,z=x+y/p,y+x/p,z-i
if d==j :
c*=math.cosh(i)
d=3*d+1
if y>=0 :
x,y,z=x-y/p,y-x/p,z+i
else :
x,y,z=x+y/p,y+x/p,z-i
j+=1
p*=2
return x*c
Inversa de un número
def inverso(n) :
a=1
r=0
for i in xrange(16) :
r+=r
if a>=n :
r+=1
a-=n
a+=a
if a>=n : r+=1
return r/32768.0
Multiplicaciones y divisiones
Consideraciones finales
