Calculadoras, senos y matrices

Después de dos entradas haciendo especial énfasis en la utilización de las matrices, creo que ya sabéis una de las patas de las que cojea este blog (pero como tenemos muchas, no nos preocupa).

Escribiendo la primera entrada acerca de las funciones generatrices me hice una pregunta muy tonta, pero que en realidad me comía por dentro…

¿cómo calculará una calculadora el seno de un ángulo?

Y pensé, cosa rara en mí, ‘…seguro que tiene implementada la serie de Taylor’.  ¡Tenía sentido! La serie de Taylor del seno en torno a x=0, con un error de no más de una centésima en el intervalo [-\pi,\pi], vimos que era,

\displaystyle f(x)=x-\frac{x^3}{3!}+\frac{x^5}{5!}-\frac{x^7}{7!}+\frac{x^9}{9!}

Aunque algo me escamaba: un algoritmo así requiere a priori un gran número de términos para obtener una buena aproximación y por tanto un buen número de multiplicaciones.

Así que empecé a buscar y encontré lo que parece ser el corazón de la mayoría de las calculadoras que se usan en la actualidad.

El algoritmo en cuestión se llama CORDICCoOrdinate Rotation DIgital Computer. Este algoritmo destaca por su simplicidad y por las operaciones necesarias para ejecutarlo: suma y bit shift, operaciones propias de la ALU de los microprocesadores. Esto hace de ORDIC un algoritmo muy rápido.

Veamos si soy capaz de explicarlo. Primero recordamos que, dado un vector en el plano u, éste puede girarse mediante la matriz,

\displaystyle R(\theta)=\begin{pmatrix}\cos\theta & -\sin\theta\\\sin\theta & \cos\theta\end{pmatrix}

Aplicando esta rotación a u, la relación entre el vector rotado, v, y u puede escribirse como,

\displaystyle \begin{cases}v_1 = \cos\theta \,u_1-\sin\theta \,u_2\\v_2 = \sin\theta \,u_1+\cos\theta \,u_2\end{cases}

que puede reescribirse también como,

\displaystyle \begin{cases}v_1 = \cos\theta \left(u_1-\tan\theta \,u_2\right)\\v_2 = \cos\theta\left(\tan\theta \,u_1+u_2\right)\end{cases}

La idea de CORDIC es aplicar sucesivas rotaciones a un vector inicial, de módulo o magnitud 1, hasta llegar al ángulo deseado. El vector inicial que suele tomarse por comodidad, y sentido común, es (1,0).

¿Por qué de módulo uno?

Porque todos los vectores de módulo 1 son de la forma,

\displaystyle (\cos\alpha,\sin\alpha)

para un cierto ángulo \alpha. Por tanto, el valor del seno, o coseno, de cualquier ángulo puede codificarse en un vector de módulo 1. 

Como comentábamos, el algoritmo consiste en realizar rotaciones. Rotaciones que no serán siempre en el mismo sentido y serán cada vez más y más pequeñas. De hecho, satisfaciendo,

\displaystyle \tan \theta_k = \pm\frac{1}{2^k}

El proceso es muy parecido al método de bisección: primero rotamos el vector inicial un ángulo de \theta_0 radianes y comparamos con el ángulo que buscamos \alpha. Si es menor, rotamos el vector en el mismo sentido que antes \theta_1 radianes, y en otro caso, rotamos en sentido contrario \theta_1 radianes…

Siguiendo con este proceso obtenemos una suma,

\displaystyle \theta_0+\theta_1+\theta_2+\theta_3-\theta_4+\theta_5-\theta_6-\theta_7-\theta_8+\theta_9+\cdots

que poco a poco (y no tan poco) se va pareciendo al valor de \alpha.

Metiendo la condición que hemos elegido para la tangente en las ecuaciones anteriores, vemos que en la rotación k+1 se tiene,

\displaystyle \begin{cases}v^{k+1}_1=\dfrac{1}{\sqrt{1+2^{-2k}}} \left(v^{k}_1-\dfrac{s}{2^k} v^{k}_2\right)\\v^{k+1}_2=\dfrac{1}{\sqrt{1+2^{-2k}}}\left(\dfrac{s}{2^k}v^{k}_1+v^{k}_2\right)\end{cases}

donde s será 1 o -1 dependiendo de si la rotación es antihoraria, o por el contrario, horaria.

¿Y cuántas rotaciones son necesarias?

Depende del nivel de detalle que queramos. Puede demostrarse que tras aplicar N+1 rotaciones se tiene,

\displaystyle \left|v^{N+1}_1 -\cos\alpha\right|\leq\frac{1}{2^N},\quad \left|v^{N+1}_2 -\sin\alpha\right|\leq \frac{1}{2^N}

Usando que,

\displaystyle 2^{-N} = 10^{-(\log_{10} 2)N} \approx 10^{-0.3 N}

vemos que tras 3 o 4 iteraciones del algoritmo se obtiene, como mínimo, una cifra significativa.

…tampoco converge muy rápido que se diga…

Toda la razón…pero como decíamos al principio, el algoritmo solo requiere de suma y bit shift (división por una potencia de 2). Aquí radica su fuerza.

Hombre…y multiplicación, ¿no?

¡No hace falta…casi! Realmente el factor \displaystyle \frac{1}{\sqrt{1+2^{-2k}}} es de normalización.

En lugar de aplicarlo en cada iteración, multiplicamos el acumulado al final del proceso. El algoritmo queda como sigue: primero las rotaciones,

\displaystyle \begin{cases}v^{k+1}_1=v^{k}_1-\dfrac{s}{2^k} v^{k}_2\\v^{k+1}_2=\dfrac{s}{2^k}v^{k}_1+v^{k}_2\end{cases}

y al final, la normalización,

\displaystyle \begin{cases}v^{N+1}_1 \leftarrow L_{N+1} v^{N+1}_1\\v^{N+1}_2 \leftarrow L_{N+1} v^{N+1}_2\end{cases}

donde \displaystyle L_{N+1}=\prod_{k=0}^{N}\frac{1}{\sqrt{1+2^{-2k}}}.

…¡entonces se multiplica!

 …¡pero sólo una vez!

Dependiendo del grado de aproximación que busquemos necesitaremos una constante de normalización u otra, pero ésta no dependerá del ángulo.

Por ejemplo, si queremos poder aproximar el seno y el coseno de un ángulo con hasta 20 cifras significativas basta guardarnos 70 constantes,

\displaystyle L_1,\ldots,L_{69}, L_{70}

y usar la adecuada en función del número de rotaciones.

Última pincelada: el algoritmo CORDIC fue desarrollado en 1959 por Jack E. Volder con fines militares y debemos comentar que su uso no se restringe al cálculo del seno y el coseno: puede adaptarse a funciones trigonométricas e hiperbólicas en general.

Espero que ahora cuando os pregunten,

…¿cómo calcula el seno una calculadora?

no digáis Taylor; decid CORDIC.

Os propongo echarle una ojeada a los siguientes artículos acerca del tema,

  1. CORDIC: How Hand Calculators Calculate de Alan Sultan.
  2. How Do Calculators Calculate Trigonometric Functions? de Jeremy M. Underwood y Bruce H. Edwards.

Esta entrada participa en El Carnaval de las Matemáticas de Abril de 2012, que tiene por anfitriones a DesEquiLIBROS. Lectura y Cultura.

Anuncios
Esta entrada fue publicada en Uncategorized y etiquetada , , , , , , , , , , , , , . Guarda el enlace permanente.

4 Responses to Calculadoras, senos y matrices

  1. serjical dice:

    Buen artículo! Chicos, no dejéis de echarle un ojo a las referencias que pone porque son magníficas. Yo en uno de mis diseños utilizo un divisor CORDIC dentro de una FPGA de Xilinx.

  2. Pingback: Calculadoras, senos y matrices

  3. DAVID dice:

    muy bueno el articulo

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s