			*Últimos Cambios*

Mi estado actual en Jabber/XMPP: - jabberES - jabber.org

El algoritmo CORDIC

Ú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.

Cálculo de senos y cosenos

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], y S[n+1]=S[n] en caso contrario. Esta opción es muy sencilla de implementar, pero tiene el problema de que Producto(cos(A[n])) no es una constante, sino que depende de los A[n] concretos que se utilicen. Ni que decir que esto es una faena.

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):

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)

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:

El número de etapas (definición de ángulos) utilizados.
Precisión de la tabla de tangentes.
Precisión de los resultados intermedios.

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.

Arcotangente y magnitud

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])

x[n+1]=sqrt(x[0]²+y[0]²)/K, donde K es el producto de los cosenos de e[n].

La rutina en Python sería:

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

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:

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

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])

x[n+1]=sqrt(x[0]²-y[0]²)/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:

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

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:

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

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.

Raíz cuadrada

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.

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

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:

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

La rutina en Python sería:

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

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.

Consideraciones finales

En muchos casos la convergencia y el error dependen del intervalo utilizado para los cálculos. En esos casos es conveniente utilizar técnicas de preescalado/postescalado.
En algunos listados propuestos aparecen multiplicaciones o divisiones como paso final antes de devolver el resultado. En muchos casos ese escalado por una constante se puede integrar en la inicialización de los vectores.
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).
Invirtiendo el movimiento de los vectores en el algoritmo para calcular logaritmos neperianos, obtendremos exponenciaciones en base e. Lo dejo como ejercicio al lector };-)