Por que o algoritmo CORDIC ficou permanentemente gravado na minha cabeça

 (github.com/francisrstokes)
3 pontos por GN⁺ 2024-05-12 | 1 comentários | Compartilhar no WhatsApp
  • CORDIC é um algoritmo que transforma cálculos complexos em operações centradas em adição e deslocamento de bits para calcular funções trigonométricas como sin, cos e tan sem FPU nem grandes tabelas de consulta
  • Esse método é mais útil em ambientes embarcados do que em sistemas de alto desempenho, especialmente em microcontroladores de baixo poder e FPGAs, e seu valor não pode ser julgado apenas pela velocidade
  • Usando ponto fixo em vez de ponto flutuante, é possível dividir os 16 bits superiores de um int32_t como parte inteira e os 16 inferiores como parte fracionária, representando aproximadamente de -32768.99997 até 32767.99997
  • Ao rotacionar um vetor com ângulos cada vez menores até o ângulo-alvo e usar 16 valores em uma tabela de atan(2**-i) junto com o valor inicial x=39796, é possível substituir a multiplicação por deslocamentos de bits em cada iteração
  • Após 16 iterações no ângulo de exemplo 0.9152, o erro absoluto de sin(0.9152) cai para 0.00000956 e o de cos(0.9152) para 0.0000434

Ambientes de cálculo em que o CORDIC se encaixa bem

  • CORDIC é um algoritmo para calcular funções trigonométricas como sin, cos e tan em hardware de baixo consumo
  • Ele funciona mesmo em ambientes sem FPU, isto é, sem unidade de ponto flutuante, ou onde é difícil usar grandes tabelas de consulta
  • O cálculo real é composto principalmente por adições e deslocamentos de bits simples
  • Ele combina matemática vetorial, trigonometria, convergência e ideias da ciência da computação para aproximar funções complexas com operações simples
  • Em hardware de alto desempenho, essa técnica pode não ser indispensável
    • Seu principal alvo de aplicação são ambientes embarcados
    • Em especial, ele é adequado para microcontroladores de baixo desempenho e FPGAs
    • Pode haver hardware ou periféricos mais rápidos, mas velocidade não é a única medida de utilidade

Representação em ponto fixo para evitar ponto flutuante

  • Mesmo funções como sin(x), que produzem valores entre -1.0 e 1.0, não precisam necessariamente ser representadas com ponto flutuante
  • Ponto fixo representa números racionais fixando a posição da vírgula dentro de um tipo inteiro
  • No exemplo, um int32_t é dividido em 16 bits superiores para a parte inteira e 16 bits inferiores para a parte fracionária
    • Nesse caso, o intervalo é aproximadamente de -32768.99997 até 32767.99997
    • Dependendo de onde a vírgula é colocada, há uma troca entre alcance da parte inteira e precisão da parte fracionária
  • O valor em si continua sendo um int32_t, e o programador atribui significado extra ao arranjo de bits

Conversão para ponto fixo e operações básicas

  • Se a precisão fracionária for de 16 bits, um valor float como 42.01 pode ser transformado em ponto fixo multiplicando por (1 << 16)
    • 42.01 * (1 << 16), ao fazer cast para int32_t, vira 2753167
    • Para voltar a float, basta calcular 2753167 / (1 << 16) e obter aproximadamente 42.0099945
  • Também é possível codificar diretamente um valor como 1.5 sem usar ponto flutuante
    • A parte inteira 1 é elevada com (1 << 16)
    • Metade da parte fracionária pode ser colocada em 0x7fff, o valor intermediário entre 0x0000 e 0xffff
    • O resultado dessa forma é o decimal 98303
  • Entre valores que usam o mesmo fator de escala, adição e subtração funcionam diretamente
  • Na multiplicação, multiplica-se os dois valores em ponto fixo e depois desloca-se o resultado para a direita pelo fator de escala
  • Na divisão, é possível obter precisão extra deslocando antes o dividendo para a esquerda pelo fator de escala e então dividindo pelo divisor

Aproximando funções trigonométricas com rotação de vetores

  • CORDIC é a sigla de “co-ordinate rotation digital computer” e foi criado em meados da década de 1950
  • A ideia central é rotacionar um vetor sobre o círculo unitário com ângulos cada vez menores, de modo que, ao atingir o ângulo-alvo, os componentes do vetor se tornem os valores de seno e cosseno
  • Esse processo avança de maneira parecida com uma busca binária
    • Move-se em direção ao ângulo-alvo com um ângulo grande
    • Verifica-se se o alvo foi ultrapassado
    • Depois repete-se a rotação no sentido horário ou anti-horário com ângulos menores
  • Como exemplo, ao calcular sin(0.7), começa-se com o vetor inicial (1, 0) e com o alvo em 0.7 radianos
    • Primeiro faz-se uma rotação anti-horária de 0.7853 radianos, isto é, 45˚
    • O alvo restante passa a ser 0.7 - 0.7853 = -0.0853
    • Como o valor é negativo, a próxima rotação é horária em 0.3926 radianos, isto é, 22.5˚
    • Depois, conforme o sinal do alvo restante, a direção vai sendo invertida com ângulos menores como 0.1963 radianos
  • Após 16 iterações, o vetor quase coincide com o ângulo original desejado, e y se torna uma aproximação de sin(a), enquanto x aproxima cos(a)

Reduzindo operações caras na matriz de rotação

  • Uma rotação de vetor comum usa multiplicação de matrizes com seno e cosseno
  • O CORDIC usa identidades trigonométricas para reescrever a matriz de rotação em torno de tan(a)
  • Como no início são usados ângulos de rotação fixos como 45˚, 22.5˚ e 11.25˚, os valores de tan(a) podem ser colocados em uma tabela pré-calculada
  • Essa tabela precisa de apenas 16 valores uint32_t, totalizando 64 bytes
    • Como comparação, uma tabela sin(x) não otimizada com 4096 valores entre -1 e 1 exigiria 16KiB e ainda teria baixa precisão
  • O termo cos(a) na frente aparece em cada rotação, mas o produto de todos eles converge para uma constante
    • Com ângulos como 45˚, 22.5˚ e 11.25˚, esse produto é aproximadamente 0.6366
    • Basta multiplicar por essa constante uma única vez ao fim de todas as iterações

Escolha dos ângulos para sobrar apenas deslocamentos e adições

  • Para eliminar multiplicações, escolhem-se ângulos cujo resultado de tan(a) seja sempre uma potência inversa de 2
  • Para isso, cria-se uma tabela com 16 entradas contendo atan(2**-i) para cada iteração i=0 até 15
  • Os ângulos reais de rotação passam a ser 45˚, 26.565˚, 14.036˚, 7.125˚ e assim por diante
  • Embora os ângulos não diminuam exatamente pela metade, o processo ainda converge corretamente usando esses valores
  • A multiplicação por tan(a) vira um deslocamento de bits pelo número de iteração i
  • O produto dos termos cos(a) também é recalculado de acordo com essa nova escolha de ângulos
    • O valor é aproximadamente 0.60725
    • Em ponto fixo de 16 bits, isso vira 39796
    • Em vez de multiplicar no final, basta iniciar o vetor com x igual a 39796 em vez de 1

Procedimento do algoritmo

  • Na etapa de pré-cálculo, monta-se uma tabela em que cada entrada é atan(2**-i) e cada valor é convertido para ponto fixo
    • A fórmula de conversão é atan(2**-i) * (1 << 16)
  • Ao calcular sin ou cos, o ângulo de entrada também é convertido para ponto fixo
    • O exemplo 0.9152 vira 0.9152 * (1 << 16) = 59978
  • O estado inicial é o seguinte
x = 39796
y = 0
z = 59978
  • z não faz parte do vetor; ele é o valor usado para acompanhar o ângulo-alvo restante
  • O sinal de z determina a direção da rotação
    • Se z >= 0, a rotação é anti-horária e executa-se z -= table[i]
    • Se z < 0, a rotação é horária e executa-se z += table[i]
  • Em cada iteração, x e y usam apenas adição, subtração e deslocamento >> i
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 convergência do exemplo e temas restantes

  • No exemplo de 0.9152 radianos, na primeira iteração z é positivo, então o vetor gira no sentido anti-horário cerca de 0.785 radianos
  • Na segunda iteração, z ainda é positivo, então ele gira cerca de 0.436 radianos no sentido anti-horário, mas passa do alvo
  • Na terceira iteração, z se torna negativo e o vetor gira cerca de 0.244 radianos no sentido horário
  • Na quarta iteração, z continua negativo, então ele gira cerca de 0.124 radianos no sentido horário
  • À medida que a variação angular diminui, o vetor oscila para frente e para trás cada vez mais perto do resultado real até convergir
  • Após 16 iterações, y se torna uma aproximação muito próxima de sin(0.9152)
    • O erro absoluto do seno é 0.00000956
    • O erro absoluto do cosseno em x é 0.0000434
  • Ainda há tópicos que não foram abordados
    • O tratamento especial necessário quando o ângulo de interesse está fora do primeiro ou do quarto quadrante do círculo unitário
    • Variações do CORDIC capazes de calcular tan, atan, asin, acos, sinh, cosh, tanh, sqrt, ln e e^x
    • O algoritmo relacionado BKM, projetado para cálculo de logaritmos e exponenciais
  • Há planos de tratar esses assuntos com mais detalhes no canal do YouTube Low Byte Productions

Leituras relacionadas

1 comentários

 
GN⁺ 2024-05-12
Opiniões do Hacker News

  • O autor disse que isso se aplica principalmente a lugares como FPGA, mas também pode ser usado em desenvolvimento de jogos ou simulações físicas distribuídas
    Cálculos de ponto flutuante tornam complicado alinhar a determinismo entre plataformas, e uma solução é evitar totalmente ponto flutuante e implementar um motor de física de ponto fixo
    Para implementar funções trigonométricas, é preciso algo como CORDIC
    Alguns anos atrás, comecei a fazer algo assim por diversão, mas não terminei, e gostaria de tentar de novo algum dia
    https://randomascii.wordpress.com/2013/07/16/floating-point-...

    • Esse texto já tem 10 anos, mas contém uma citação importante: “O padrão IEEE garante algumas coisas. Garante mais do que sabem as pessoas que veem a matemática de ponto flutuante como uma espécie de misticismo, mas menos do que alguns programadores imaginam”
      Em resumo, o x87 tinha suas esquisitices; é preciso manter consistentes configurações como modo de arredondamento e flush-to-zero; processadores antigos não têm FMA; instruções aproximadas como mmsqrtps não têm uma especificação consistente; e o compilador pode reassociar expressões
      Em uma rotina pequena ou em uma biblioteca escrita por você, embora seja doloroso, é possível garantir que essas coisas sejam evitadas
      O IEEE-754 2008 deixou a especificação mais clara e, na prática, partiu do pressuposto da morte do x87; em 2024, dá para evitar o x87 com segurança
      FMA também faz parte da especificação IEEE-754 2008 e está presente em processadores modernos, incluindo os posteriores ao Intel Haswell
      Ainda assim, diferenças de arquitetura como AVX2 de 8 lanes e NEON de 4 lanes podem atrapalhar, mas, usando assembly, intrínsecos ou C verificado com Compiler Explorer ou objdump, dá para olhar a saída e julgar: “isso será consistente”
    • O próprio autor também disse que, antes de o ponto flutuante se popularizar graças ao aumento de desempenho do hardware, ponto fixo era muito comum no desenvolvimento de jogos, e é bem provável que CORDIC também fosse usado junto
      “Na verdade, antes de o IEEE 754 se tornar um padrão tão popular quanto é hoje, ponto fixo era usado o tempo todo. Pergunte a um desenvolvedor de jogos que trabalhou entre 1980 e por volta de 2000, e ele contará os detalhes”
    • A biblioteca de simulação física para desenvolvimento de jogos nphysics adotou uma abordagem que usa matemática de ponto fixo junto com CORDIC quando é necessário determinismo entre plataformas, mas hoje foi descontinuada
      A nova biblioteca Rapier, reescrita a partir da nphysics, em vez disso se apoia nas garantias do IEEE-754 2008 para oferecer determinismo entre plataformas
      Por isso, não funciona em plataformas antigas, mas é determinística em plataformas modernas, incluindo wasm
      É claro que não dá para depender das rotinas de funções transcendentais como sin e cos fornecidas por cada plataforma; é preciso implementá-las diretamente para que se comportem da mesma forma em todos os lugares
      Mas, se não for executada em plataformas não conformes, é uma abordagem viável
      https://www.rustsim.org/blog/2020/06/01/this-month-in-rustsi...
      https://rapier.rs/docs/user_guides/rust/determinism/

  • CORDIC pode ser usado não só para calcular e gerar seno e cosseno, mas também para várias operações como logaritmos, exponenciais, raízes quadradas, magnitude de vetores, conversão entre coordenadas polares e cartesianas e rotação de vetores
    O autor também antecipa essas possibilidades na conclusão
    Tenho a impressão de que, usando quatérnios em vez das matrizes ortonormais convencionais, operações baseadas em CORDIC poderiam ser executadas de forma mais eficiente, isto é, usando menos ciclos de computação e menos memória, além de reduzir erros
    https://core.ac.uk/works/8439118

    • Se bem me lembro, também pode ser estendido para grupos de Lie arbitrários

  • Aprendi séries de Taylor no pré-cálculo do ensino médio, e o professor disse que as funções trigonométricas das calculadoras eram implementadas assim na prática
    Quando fui pesquisar, vi que na verdade era CORDIC, e me diverti implementando em TI Basic

    • Talvez seja interessante ler sobre como a incrível calculadora Sinclair Scientific calculava funções trigonométricas, logaritmos e afins
      Não era CORDIC, mas há semelhanças no algoritmo
      http://files.righto.com/calculator/sinclair_scientific_simul...
    • Existe mesmo alguma calculadora que use expansão de Taylor?

  • Textos sobre implementação em hardware:
    https://arxiv.org/pdf/2211.04053
    https://hal.science/hal-01327460/document
    https://archive.ll.mit.edu/HPEC/agendas/proc05/Day_1/Abstrac...
    Gostaria de ver como isso se compara, em diferentes hardwares ao longo do tempo, às implementações de funções trigonométricas comuns em software e hardware

    • É estranho que o CORDIC, embora seja uma técnica computacional amplamente usada e muito popular, não seja tratado em detalhes adequadamente nos livros
      Com o crescimento de IoT e da comunicação máquina a máquina, e considerando as implementações de CORDIC e sua eficiência computacional, seu uso provavelmente deve aumentar bastante; portanto, são necessárias boas referências para implementações corretas e otimizadas
      Como exceção, há os livros do Prof. Omondi e do Prof. Deschamps
      https://www.worldscientific.com/worldscibooks/10.1142/p1054
      http://www.arithmetic-circuits.org/guide2fpga/vhdl_codes.htm

  • sin e cos são usados com frequência em rotações de vetores
    Nesse caso, o truque do CORDIC é evitar o cálculo tradicional de sin/cos/multiplicação e passar o próprio vetor a ser rotacionado como entrada do CORDIC
    Assim, o CORDIC produz diretamente o vetor rotacionado sem calcular sin/cos nem fazer multiplicação de números complexos
    O CORDIC brilha especialmente quando a latência não é tão importante
    Se cada etapa do cálculo for colocada em pipeline, é possível obter uma vazão alta, o que combina bem com mixagem digital em sistemas sem fio

  • Em 2023, alguns MCUs modernos são baratos e ainda assim têm FPU
    O STM32G4 é um bom exemplo e, diferentemente de casos como MCUs M0, se você não quiser usar ponto fixo, pode usar f32 à vontade
    Esses chips podem ser encontrados por cerca de 1 a 2 dólares por MCU
    Ainda assim, o G4 também tem um periférico CORDIC em hardware que implementa esse algoritmo para uso com ponto fixo
    Fico curioso se isso serve principalmente para evitar perdas de precisão de ponto flutuante
    Ele é programado por registradores, mas não é uma implementação direta do CORDIC na CPU; um hardware dedicado dentro do IC é que faz o processamento

  • Isso me lembrou um trecho de código bem bonitinho em que participei no passado
    Era preciso encontrar as coordenadas da bissetriz do ângulo formado por um arco no círculo unitário, e as coordenadas (x,y) dos dois braços já estavam disponíveis
    A implementação existente era um amontoado de trigonometria: convertia as coordenadas (x,y) para coordenadas polares (r,θ), verificava se o θ calculado estava no quadrante correto, dividia θ pela metade e então convertia de volta para (x,y)
    No fim, chamava muitas funções trigonométricas e inversas
    Como era Python e números complexos podiam ser usados como cidadãos de primeira classe, bastou definir dois números complexos, z1 a partir de (x1,y1) e z2 a partir de (x2,y2), e tirar apenas a média geométrica do produto √(z1*z2)
    No código novo não havia funções trigonométricas explícitas, nem conversões e conversões inversas explícitas

  • Está escrito: “é bastante óbvio que girar 22,75˚ é o mesmo que girar 45˚ e depois -22,5˚”; então não seria uma rotação de 22,5°?
    Fico na dúvida se é um erro do texto ou se eu entendi errado

    • É um erro do texto

  • O sistema de octree de Meagher é famoso por usar apenas aritmética inteira, sem multiplicação nem divisão de inteiros
    “Foram desenvolvidos algoritmos eficientes de tempo linear para exibição, incluindo operações booleanas (união, interseção, diferença), operações geométricas (translação, redimensionamento, rotação), detecção de interferência em N dimensões e remoção de superfícies ocultas em pontos arbitrários do espaço. Esses algoritmos não requerem operações de ponto flutuante, multiplicação de inteiros nem divisão de inteiros”
    https://doi.org/10.1016/0146-664X(82)90104-6
    Graças a isso, ficou mais fácil criar hardware acelerador gráfico VLSI rápido e personalizado para representações em octree

  • Fico curioso sobre o desempenho do CORDIC em comparação com interpolação cúbica usando uma tabela pequena, ou com outras interpolações polinomiais
    Aprendi que sintetizadores com recursos limitados às vezes usavam interpolação cúbica; provavelmente era uma época em que o CORDIC ainda era relativamente novo
    Olhando por alto, como o CORDIC ganha 1 bit de precisão a cada iteração, o cálculo deve ser mais caro, mas deve usar menos espaço do que polinômios
    Porém, do ponto de vista de espaço, é importante destacar que ele pode sair mais barato do que a tabela de consulta com 4096 itens para sin(x) apresentada no texto
    Graças à simetria, só é necessário 1/4 do círculo inteiro

    • Antigos desenvolvedores de jogos e da demoscene usavam tabelas de consulta de apenas 256 itens para sin e cos
      Usar ângulos do tamanho de um byte era conveniente porque eles circulavam automaticamente, e 2^8 era bastante suficiente para rotações em jogos 2D
      Mas, se você quiser movimentos suaves, em 3D isso não vai muito longe