Por qué el algoritmo CORDIC quedó grabado permanentemente en mi cabeza

 (github.com/francisrstokes)
3 puntos por GN⁺ 2024-05-12 | 1 comentarios | Compartir por WhatsApp
  • CORDIC es un algoritmo que, para calcular funciones trigonométricas como sin, cos y tan sin una FPU ni grandes tablas de consulta, reemplaza operaciones complejas por un enfoque centrado en sumas y desplazamientos de bits
  • Este método es útil no tanto en sistemas de alto rendimiento como en entornos embebidos, especialmente en microcontroladores de bajo desempeño y FPGA, donde es difícil juzgar su valor solo por la velocidad
  • Si se usa punto fijo en lugar de punto flotante, se pueden dividir los 32 bits de int32_t en 16 bits superiores para la parte entera y 16 bits inferiores para la parte fraccionaria, con lo que se puede representar aproximadamente desde -32768.99997 hasta 32767.99997
  • Al rotar un vector hacia un ángulo objetivo con incrementos cada vez más pequeños, y usando una tabla de 16 valores atan(2**-i) junto con un valor inicial x=39796, se puede sustituir la multiplicación de cada iteración por desplazamientos de bits
  • Si se repite 16 veces con el ángulo de ejemplo 0.9152, el error absoluto de sin(0.9152) se reduce a 0.00000956 y el de cos(0.9152) a alrededor de 0.0000434

Entorno de cálculo ideal para CORDIC

  • CORDIC es un algoritmo para calcular funciones trigonométricas como sin, cos y tan en hardware de bajo consumo
  • Funciona incluso en entornos donde no hay FPU, es decir, una unidad de punto flotante, o donde es difícil usar tablas de consulta grandes
  • El cálculo real se compone principalmente de sumas y desplazamientos de bits
  • Combina matemáticas vectoriales, trigonometría, convergencia e ideas de ciencias de la computación para aproximar funciones complejas con operaciones simples
  • En hardware de alto rendimiento, esta técnica puede no ser indispensable
    • Su principal ámbito de aplicación es el entorno embebido
    • En particular, es adecuada para microcontroladores de bajo desempeño y FPGA
    • Puede haber hardware o periféricos más rápidos, pero la velocidad no es la única medida de utilidad

Representación en punto fijo para evitar el punto flotante

  • Incluso una función como sin(x), que produce valores entre -1.0 y 1.0, no necesita expresarse obligatoriamente en punto flotante
  • El punto fijo representa números racionales fijando la posición del punto decimal dentro de un tipo entero
  • El ejemplo divide int32_t en 16 bits superiores para la parte entera y 16 bits inferiores para la parte fraccionaria
    • En ese caso, el rango va aproximadamente de -32768.99997 a 32767.99997
    • Según dónde se coloque el punto decimal, se intercambia alcance en la parte entera por precisión en la parte fraccionaria
  • El valor en sí sigue siendo un int32_t, y es el programador quien le da significado adicional al arreglo de bits

Conversión a punto fijo y operaciones básicas

  • Si la precisión fraccionaria es de 16 bits, un valor float como 42.01 puede convertirse a punto fijo multiplicándolo por (1 << 16)
    • 42.01 * (1 << 16) se convierte en 2753167 al hacer cast a int32_t
    • Para volver a float, basta calcular 2753167 / (1 << 16) y se obtiene aproximadamente 42.0099945
  • También se puede codificar directamente un valor como 1.5 sin usar punto flotante en absoluto
    • La parte entera 1 se sube con (1 << 16)
    • La mitad fraccionaria puede colocarse como 0x7fff, el valor medio entre 0x0000 y 0xffff
    • El resultado de este método es el decimal 98303
  • Entre valores que usan el mismo factor de escala, la suma y la resta funcionan sin cambios
  • La multiplicación se hace multiplicando dos valores en punto fijo y luego desplazando el resultado a la derecha por el factor de escala
  • En la división, si primero se desplaza el dividendo a la izquierda por el factor de escala y luego se divide por el divisor, se obtiene precisión adicional

Aproximar funciones trigonométricas con rotación de vectores

  • CORDIC significa “co-ordinate rotation digital computer” y fue creado a mediados de los años 50
  • La idea central es rotar un vector sobre el círculo unitario con ángulos cada vez más pequeños, de modo que al llegar al ángulo objetivo los componentes del vector sean los valores seno y coseno
  • Este proceso avanza de forma parecida a una búsqueda binaria
    • Se avanza con un ángulo grande hacia el ángulo objetivo
    • Se comprueba si ya se pasó del objetivo
    • Luego se repiten rotaciones en sentido horario o antihorario con ángulos más pequeños
  • Por ejemplo, para calcular sin(0.7), se empieza con el vector inicial (1, 0) y el objetivo de 0.7 radianes
    • Primero se rota 0.7853 radianes, es decir, 45˚, en sentido antihorario
    • El objetivo restante pasa a ser 0.7 - 0.7853 = -0.0853
    • Como el valor es negativo, la siguiente rotación es de 0.3926 radianes, es decir, 22.5˚, en sentido horario
    • Después, según el signo del objetivo restante, se sigue rotando con ángulos más pequeños como 0.1963 radianes, cambiando de dirección
  • Tras 16 iteraciones, el vector queda casi alineado con el ángulo objetivo original, y y es una aproximación de sin(a) mientras que x es una aproximación de cos(a)

Reducir operaciones costosas en la matriz de rotación

  • La rotación de vectores convencional usa multiplicación de matrices con seno y coseno
  • CORDIC transforma la matriz de rotación para centrarla en tan(a) usando identidades trigonométricas
  • Al principio se usan ángulos de rotación fijos como 45˚, 22.5˚ y 11.25˚, por lo que los valores de tan(a) pueden precomputarse en una tabla
  • Esta tabla solo necesita 16 valores uint32_t, es decir, 64 bytes
    • Como comparación, una tabla sin(x) no optimizada con 4096 valores entre -1 y 1 requiere 16KiB y además se considera menos precisa
  • El término cos(a) que aparece al frente en cada rotación se genera en cada iteración, pero el producto de todos esos términos converge a una constante
    • Para ángulos como 45˚, 22.5˚ y 11.25˚, ese producto es aproximadamente 0.6366
    • Esa constante puede multiplicarse una sola vez al final de todas las iteraciones

Elegir ángulos para quedarse solo con desplazamientos y sumas

  • Para eliminar la multiplicación, se eligen ángulos tales que el resultado de tan(a) sea siempre una potencia inversa de 2
  • Para ello, se crea una tabla de 16 entradas con los valores atan(2**-i) para cada iteración i=0 a 15
  • Los ángulos de rotación reales pasan a ser 45˚, 26.565˚, 14.036˚, 7.125˚, etc.
  • Los ángulos no se reducen exactamente a la mitad cada vez, pero aun así el proceso converge al resultado correcto
  • La multiplicación por tan(a) se convierte en un desplazamiento de bits según el número de iteración i
  • El producto de los términos cos(a) también se recalcula de acuerdo con esta nueva elección de ángulos
    • Su valor es aproximadamente 0.60725
    • En punto fijo de 16 bits, eso equivale a 39796
    • En vez de multiplicar al final, basta con fijar la x inicial en 39796 en lugar de 1

Procedimiento del algoritmo

  • En la etapa de precálculo se construye una tabla cuyas entradas son atan(2**-i) y luego cada valor se convierte a punto fijo
    • La fórmula de conversión es atan(2**-i) * (1 << 16)
  • Para obtener sin o cos, el ángulo de entrada también se convierte a punto fijo
    • El ejemplo 0.9152 se convierte en 0.9152 * (1 << 16) = 59978
  • El estado inicial es el siguiente
x = 39796
y = 0
z = 59978
  • z no es parte del vector, sino un valor que sigue el ángulo objetivo restante
  • El signo de z determina la dirección de rotación
    • Si z >= 0, se rota en sentido antihorario y se hace z -= table[i]
    • Si z < 0, se rota en sentido horario y se hace z += table[i]
  • Cada iteración usa solo suma, resta y desplazamientos >> i sobre x y y
if z >= 0:
    x_next = x - (y >> i)
    y_next = y + (x >> i)
    z -= table[i]
else:
    x_next = x + (y >> i)
    y_next = y - (x >> i)
    z += table[i]
x = x_next
y = y_next

Resultado de convergencia del ejemplo y temas pendientes

  • En el ejemplo de 0.9152 radianes, en la primera iteración z es positivo, así que se rota en sentido antihorario unos 0.785 radianes
  • En la segunda iteración z sigue siendo positivo, así que se rota otros 0.436 radianes en sentido antihorario, pero se sobrepasa el objetivo
  • En la tercera iteración z se vuelve negativo, así que se rota en sentido horario unos 0.244 radianes
  • En la cuarta iteración z también es negativo, así que se rota en sentido horario unos 0.124 radianes
  • A medida que el cambio angular se hace más pequeño, el vector oscila hacia adelante y atrás cerca del resultado real hasta converger
  • Después de 16 iteraciones, y se vuelve una aproximación muy cercana de sin(0.9152)
    • El error absoluto del seno es 0.00000956
    • El error absoluto del coseno en x es 0.0000434
  • Quedan temas que no se cubrieron
    • El tratamiento especial necesario cuando el ángulo de interés queda fuera del primer o cuarto cuadrante del círculo unitario
    • Otras funciones que pueden calcularse con variantes de CORDIC: tan, atan, asin, acos, sinh, cosh, tanh, sqrt, ln, e^x
    • El algoritmo relacionado BKM, diseñado para cálculos de logaritmos y exponenciales
  • Está previsto tratar estos temas con más detalle en el canal de YouTube de Low Byte Productions

Lecturas relacionadas

1 comentarios

 
GN⁺ 2024-05-12
Comentarios de Hacker News

  • El autor dijo que se aplica sobre todo en cosas como FPGA, pero también puede usarse en desarrollo de videojuegos o simulaciones físicas distribuidas.
    Hacer que los cálculos de punto flotante sean deterministas entre plataformas es complicado, y una solución es evitar por completo el punto flotante e implementar un motor de física de punto fijo.
    Para implementar funciones trigonométricas, se necesita algo como CORDIC.
    Hace unos años empecé a hacer algo así por diversión, pero no lo terminé, y algún día me gustaría retomarlo.
    https://randomascii.wordpress.com/2013/07/16/floating-point-...

    • Ese artículo ya tiene 10 años, pero incluye una cita importante: “El estándar IEEE garantiza algunas cosas. Garantiza más de lo que saben quienes ven las matemáticas de punto flotante como una especie de misticismo, pero garantiza menos de lo que creen algunos programadores”.
      En resumen, x87 tenía rarezas, hay que mantener consistentes configuraciones como el modo de redondeo y flush-to-zero, los procesadores antiguos no tenían FMA, instrucciones aproximadas como mmsqrtps no tienen una especificación consistente, y el compilador puede reassociar expresiones.
      En rutinas pequeñas o bibliotecas escritas a mano, aunque sea doloroso, es posible asegurarse de evitar estas cosas.
      IEEE-754 2008 aclaró más la especificación y de hecho asumió la muerte de x87, y en 2024 se puede evitar x87 con seguridad.
      FMA también es parte de la especificación de IEEE-754 2008 y está presente en procesadores modernos, incluido Intel desde Haswell.
      Aun así, diferencias de arquitectura como 8-wide AVX2 y 4-wide NEON pueden causar problemas, pero usando ensamblador, intrinsics, o C verificado con Compiler Explorer u objdump, se puede mirar la salida y concluir: “esto debería ser consistente”.
    • El autor también dijo que antes de que el punto flotante se volviera ampliamente usado en desarrollo de hardware por las mejoras de rendimiento, el punto fijo era muy común en desarrollo de videojuegos, y es muy probable que CORDIC también se usara junto con ello.
      “De hecho, antes de que IEEE 754 se convirtiera en un estándar tan popular como lo es hoy, siempre se usaba punto fijo. Pregúntale a cualquier desarrollador de juegos que haya trabajado entre 1980 y más o menos 2000 y te lo contará en detalle”.
    • La biblioteca de simulación física para desarrollo de videojuegos nphysics eligió usar matemáticas de punto fijo junto con CORDIC cuando se necesitaba determinismo entre plataformas, pero ahora está descontinuada.
      Rapier, la nueva biblioteca que reescribió nphysics, en cambio ofrece determinismo entre plataformas apoyándose en las garantías de IEEE-754 2008.
      Por eso no funciona en plataformas antiguas, pero sí es determinista en plataformas modernas, incluido wasm.
      Claro, no se puede depender de rutinas de funciones trascendentales como sin o cos que provee cada plataforma, y hay que implementarlas directamente para que funcionen igual en todas partes.
      Pero es un enfoque viable si no se ejecuta en plataformas no conformes.
      https://www.rustsim.org/blog/2020/06/01/this-month-in-rustsi...
      https://rapier.rs/docs/user_guides/rust/determinism/

  • CORDIC puede usarse no solo para calcular y generar seno y coseno, sino también para varias operaciones como logaritmos, exponentes, raíces cuadradas, magnitud de vectores, conversión entre coordenadas polares y cartesianas, y rotación de vectores.
    El autor también adelanta estas posibilidades en la conclusión.
    Da la impresión de que, si se usan cuaterniones en lugar de matrices ortonormales convencionales, las operaciones basadas en CORDIC podrían ejecutarse de forma más eficiente, es decir, usando menos ciclos de cómputo y memoria, y reduciendo además el error.
    https://core.ac.uk/works/8439118

    • Si no recuerdo mal, también puede extenderse a grupos de Lie arbitrarios.

  • En precálculo de preparatoria aprendí sobre las series de Taylor, y el profesor dijo que así era como realmente se implementaban las funciones trigonométricas en las calculadoras.
    Luego investigué y en realidad era CORDIC, y me divertí bastante implementándolo en TI Basic.

    • Probablemente sería interesante leer cómo la asombrosa calculadora Sinclair Scientific calculaba funciones trigonométricas y logaritmos, entre otras cosas.
      No era CORDIC, pero el algoritmo tiene ciertas similitudes.
      http://files.righto.com/calculator/sinclair_scientific_simul...
    • ¿Existe de verdad хотя sea una sola calculadora que use desarrollos de Taylor?

  • Artículos sobre implementación en hardware:
    https://arxiv.org/pdf/2211.04053
    https://hal.science/hal-01327460/document
    https://archive.ll.mit.edu/HPEC/agendas/proc05/Day_1/Abstrac...
    Quisiera ver cómo se compara con implementaciones generales de funciones trigonométricas en software y hardware en distintos tipos de hardware a lo largo del tiempo

    • CORDIC es una técnica informática muy usada y muy popular, así que resulta extraño que los libros no la traten realmente en detalle
      Como el IoT y la comunicación máquina a máquina están creciendo, y considerando la implementación de CORDIC y su eficiencia computacional, probablemente su uso aumente mucho, por lo que hacen falta buenas referencias para una implementación correcta y optimizada
      Como excepción, están los libros del Prof. Omondi y el Prof. Deschamps
      https://www.worldscientific.com/worldscibooks/10.1142/p1054
      http://www.arithmetic-circuits.org/guide2fpga/vhdl_codes.htm

  • sin y cos se usan con frecuencia para la rotación de vectores
    En este caso, el truco de CORDIC es evitar el cálculo tradicional de sin/cos/multiplicación y meter como entrada a CORDIC el propio vector que se quiere rotar
    Entonces CORDIC produce directamente el vector rotado sin calcular sin/cos ni hacer multiplicación compleja
    CORDIC brilla especialmente cuando la latencia no es demasiado importante
    Si se canaliza cada etapa del cálculo, se puede obtener un gran rendimiento, así que encaja bien para mezcla digital en sistemas inalámbricos

  • A partir de 2023, algunos MCU modernos tienen FPU aunque sean de bajo costo
    El STM32G4 es un buen ejemplo y, a diferencia de casos como los MCU M0, si no quieres usar punto fijo puedes usar f32 sin problema
    Estos chips pueden conseguirse por alrededor de 1 a 2 dólares por MCU
    Sin embargo, el G4 también tiene un periférico CORDIC por hardware que implementa este algoritmo para punto fijo
    Me pregunto si esto se usa principalmente para evitar la pérdida de precisión del punto flotante
    Se programa mediante registros, pero no es una implementación directa de CORDIC en la CPU, sino que lo maneja hardware dedicado dentro del IC

    • Según el inventario de Digi-Key, sin contar duplicados, los Cortex-M4F más baratos son el Nuvoton M481LE8AE de 3 dólares https://www.digikey.com/en/products/detail/nuvoton-technolog..., el Maxim MAX32660 de 3 dólares https://www.digikey.com/en/products/detail/analog-devices-in... y el Atmel ATSAMD51 de 5 dólares https://www.digikey.com/en/products/detail/microchip-technol...
      El STM32G4 más barato es el STM32G441KBT6 y redondeando cuesta 4 dólares https://www.digikey.com/en/products/detail/microchip-technol...
      Me pregunto dónde se consiguen por menos de 2 dólares
      En Digi-Key, los chips de Nuvoton apenas quedan por debajo de 2 dólares en cantidades de 500 unidades
    • El segundo chip Parallax Propeller tiene un motor CORDIC implementado en silicio
      Es rápido y maneja productos intermedios de 64 bits, así que la precisión de la división y de las funciones trigonométricas es suficiente para la mayoría de los usos
      Si hace falta, la precisión puede mejorarse más por software
      Conocí CORDIC bastante tarde; antes de eso usaba mucho punto fijo en el mundo del ensamblador de 8 y 16 bits por rendimiento y determinismo
      Cuando lo conocí, me sorprendió
      Era rápido, y las matemáticas necesarias para usarlo de forma útil eran solo las básicas

  • Esto me recuerda un fragmento de código bastante simpático en el que participé hace tiempo
    Había que encontrar las coordenadas de la bisectriz del ángulo formado por un arco de la circunferencia unitaria, y ya se tenían las coordenadas (x,y) de ambos brazos
    La implementación existente era una maraña trigonométrica que convertía las coordenadas (x,y) a coordenadas polares (r,θ), comprobaba que el θ calculado estuviera en el cuadrante correcto, luego dividía θ entre dos y volvía a convertir a (x,y)
    Al final, hacía muchas llamadas a funciones trigonométricas y a sus inversas
    Como en Python los números complejos pueden usarse como ciudadanos de primera clase, bastaba con definir dos complejos z1 a partir de (x1,y1) y z2 a partir de (x2,y2), y luego tomar la media geométrica del producto √(z1*z2)
    El código nuevo no tenía trigonometría explícita, ni conversiones explícitas ni conversiones inversas

  • Dice: “rotar 22.75˚ es bastante obviamente lo mismo que rotar 45˚ y luego -22.5˚”, pero entonces, ¿no sería una rotación de 22.5°?
    Me pregunto si es un error del texto o si entendí mal

    • Es un error del texto

  • El sistema de octree de Meagher era conocido por usar solo aritmética entera, sin multiplicación ni división entera
    “Se desarrollaron algoritmos eficientes de tiempo lineal para operaciones booleanas (unión, intersección, diferencia), operaciones geométricas (traslación, cambio de escala, rotación), detección de interferencias en N dimensiones y visualización, incluida la eliminación de superficies ocultas en cualquier punto arbitrario del espacio. Estos algoritmos no requieren operaciones de punto flotante, multiplicación entera ni división entera”
    https://doi.org/10.1016/0146-664X(82)90104-6
    Gracias a esto, era fácil crear hardware acelerador gráfico VLSI rápido y personalizado para la representación con octrees

  • Me da curiosidad qué rendimiento tiene CORDIC frente a la interpolación cúbica que usa tablas pequeñas u otras interpolaciones polinomiales
    Aprendí que los sintetizadores con recursos limitados a veces usaban interpolación cúbica, aunque probablemente era una época en la que CORDIC todavía era relativamente nuevo
    A grandes rasgos, CORDIC gana 1 bit de precisión por iteración, así que el cálculo parece más costoso, pero probablemente use menos espacio que los polinomios
    Aun así, en términos de espacio hay que recalcar que puede salir más barato de lo que sugiere en el artículo la tabla de consulta de 4096 entradas para sin(x)
    Gracias a la simetría, solo hace falta 1/4 del círculo completo

    • Los antiguos desarrolladores de videojuegos y de la demoscene usaban tablas de consulta de apenas 256 entradas para sin y cos
      Si usabas ángulos de un byte, el ciclo era automático y práctico, y para la rotación en juegos 2D 2^8 era bastante suficiente
      Pero si querías movimiento suave, en 3D eso no llegaba muy lejos