Trigonometría de punto fijo para aplicaciones integradas

Question

Necesito hacer transformaciones rotacionales (y otras) en una aplicación embebida, requiriendo las funciones sin() cos() y tan(). Sé que se pueden utilizar tablas de búsqueda, y esa es la única solución que he podido encontrar investigando por mi cuenta, pero ¿hay alguna buena biblioteca de trigonometría de punto fijo por ahí?

Estoy pensando en utilizar un cortex M3 para la aplicación, por lo que quiero mantenerme alejado del punto flotante tanto como sea posible para mantener las aplicaciones zippy.

Answer 1

Un buen enfoque para realizar la trigonometría en aplicaciones embebidas es utilizar aproximaciones polinómicas a las funciones que necesitas. El código es compacto, los datos constan de unos pocos coeficientes y las únicas operaciones necesarias son multiplicar y sumar/restar. Muchos sistemas embebidos disponen de multiplicadores por hardware, que ofrecen un buen rendimiento.

Answer 2

¿Te opones a utilizar las bibliotecas Cortex de punto fijo para esto?

q31_t arm_sin_q31 (q31_t x)
Aproximación rápida a la función seno trigonométrica para datos Q31.

de:

CMSIS-DSP: Colección de librerías DSP con más de 60 funciones para varios tipos de datos: punto fijo (fraccionario q7, q15, q31) y punto flotante de precisión simple (32 bits). La biblioteca está disponible para Cortex-M0, Cortex-M3 y Cortex-M4.

Utiliza una tabla de búsqueda con interpolación cuadrática, pero es bastante rápido. Podrías adaptarlo a la interpolación lineal para obtener una mayor velocidad pero más errores.

También hay que tener en cuenta que incluso el Cortex M4 no tiene necesariamente FPU. He visto que los llaman "M4F" si lo tienen.

Answer 3

Esta respuesta pretende aumentar la respuesta actualmente aceptada con un ejemplo concreto en dos variantes, y proporcionar algunos consejos de diseño específicos.

La aproximación polinómica suele ser un enfoque superior si la precisión deseada es bastante alta y se dispone de un multiplicador de hardware. El tamaño de las tablas tiende a aumentar rápidamente incluso cuando se utilizan esquemas de interpolación (por ejemplo, lineal, cuadrática) y compresión (por ejemplo, tablas bipartitas) una vez que se requieren más de unos 16 bits "buenos".

El uso de aproximaciones minimax para los polinomios es muy recomendable, ya que minimizan el error máximo en el intervalo para el que se generan. Esto puede conducir a una reducción significativa en el número de términos requeridos para una precisión particular en comparación con, por ejemplo, las expansiones de la serie de Taylor que proporcionan la mejor precisión sólo en el punto alrededor del cual se expanden. Herramientas de uso común como Mathematica, Maple y el programa de código abierto Herramienta Sollya ofrecen métodos incorporados para generar aproximaciones minimax.

Las operaciones de multiplicación por lo alto son el bloque computacional fundamental de la evaluación de polinomios en la aritmética de punto fijo. Devuelven la mitad más significativa del producto completo de una multiplicación de enteros. La mayoría de las arquitecturas proporcionan variantes con y sin signo, otras proporcionan multiplicaciones con resultados de doble ancho devueltos en dos registros. Algunas arquitecturas ofrecen incluso combinaciones de multiplicación-agregación, que pueden ser especialmente útiles. Los compiladores de optimización suelen ser capaces de traducir los modismos del código fuente HLL (como los utilizados en el código ISO-C que se muestra a continuación) correspondientes a estas operaciones a las instrucciones hardware adecuadas.

Para maximizar la precisión de la evaluación de polinomios, se querría utilizar el máximo número de bits posible en todo momento durante el cálculo intermedio, eligiendo un formato de punto fijo con el máximo número posible de bits de fracción. Para la eficiencia, un factor de escala igual al ancho del registro evita la necesidad de re-escalar a través de desplazamientos cuando se utiliza junto con las operaciones de multiplicación-alta.

Mientras que el esquema de Horner se utiliza normalmente en el cálculo de punto flotante para evaluar polinomios de alta precisión, esto es a menudo innecesario en el cálculo de punto fijo y puede ser perjudicial para el rendimiento debido a la larga cadena de dependencia de la evaluación de polinomios que expone la latencia de la multiplicación. Suelen ser preferibles los esquemas de evaluación en paralelo que permiten la mejor utilización de los multiplicadores canalizados con latencia de varios ciclos. En el código de abajo combino los términos de cada polinomio por pares y construyo la evaluación del polinomio completo a partir de eso.

El código ISO-C que se muestra a continuación demuestra el cálculo simultáneo del seno y el coseno según estos principios de diseño en el intervalo [0, /2], donde la entrada y la salida están en formato S8.23 (Q8.23). Se obtienen resultados totalmente precisos, con un error máximo del orden del 10 -7 y el 80+% de los resultados se redondean correctamente.

La primera variante, en sincos_fixed_nj() utiliza un enfoque clásico de reducción del argumento en [0, /4], y la aproximación polinómica al seno y al coseno en ese intervalo. A continuación, la fase de reconstrucción asigna los valores polinómicos al seno y al coseno en función del cuadrante. La segunda variante, en sincos_fixed_ollyw se basa en un entrada del blog de OllyW . Proponen aplicar la transformada a = (2/)x-1/2 en el intervalo [-1/2, 1/2], en el que luego hay que aproximar el sen ((2a + )/4 y el cos ((2a + )/4. Las expansiones en serie de estos ( sin , porque ) son idénticos, salvo que el signo se invierte para los términos de potencia impar. Esto significa que uno puede sumar los términos de potencia impar y potencia par por separado y luego calcular el seno y el coseno como la suma y la diferencia de las sumas acumuladas.

Utilizando Explorador de compiladores He compilado con Clang 11.0 para un armv7-a Objetivo ARM de 32 bits con optimización completa ( -O3 ). Ambas variantes se compilan en 41 instrucciones subrutinas con cada subrutina utilizando nueve constantes de 32 bits almacenadas. sincos_fixed_ollyw() utiliza una instrucción de multiplicación más que sincos_fixed_nj pero tiene una presión de registro ligeramente inferior. La situación parece ser similar cuando se construye con Clang para otros objetivos de arquitectura, por lo que uno querría probar ambas variantes para ver cuál funciona mejor en una plataforma determinada. La tangente podría calcularse dividiendo el resultado del seno por el resultado del coseno.

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <math.h>

#define SINCOS_NJ    (1)
#define SINCOS_OLLYW (2)
#define VARIANT      (SINCOS_NJ)

/* a single instruction in many 32-bit architectures */
uint32_t umul32hi (uint32_t a, uint32_t b)
{
    return (uint32_t)(((uint64_t)a * b) >> 32);
}

/* a single instruction in many 32-bit architectures */
int32_t mul32hi (int32_t a, int32_t b)
{
    return (int32_t)(uint32_t)((uint64_t)((int64_t)a * b) >> 32);
}

/*
  compute sine and cosine of argument in [0, PI/2]
  input and output in S8.23 format
  max err sine = 9.86237533e-8  max err cosine = 1.02729891e-7
  rms err sine = 4.11141973e-8  rms err cosine = 4.11752018e-8
  sin correctly rounded: 10961278 (83.19%)  
  cos correctly rounded: 11070113 (84.01%)
*/
void sincos_fixed_nj (int32_t x, int32_t *sine, int32_t *cosine)
{
    // minimax polynomial approximation for sine on [0, PI/4]
    const uint32_t s0 = (uint32_t)(1.9510998390614986e-4 * (1LL << 32) + 0.5);
    const uint32_t s1 = (uint32_t)(8.3322080317884684e-3 * (1LL << 32) + 0.5);
    const uint32_t s2 = (uint32_t)(1.6666648373939097e-1 * (1LL << 32) + 0.5);
    const uint32_t s3 = (uint32_t)(9.9999991734512150e-1 * (1LL << 32) + 0.5);
    // minimax polynomial approximation for cosine on [0, PI/4]
    const uint32_t c0 = (uint32_t)(1.3578890357166529e-3 * (1LL << 32) + 0.5);
    const uint32_t c1 = (uint32_t)(4.1654359549283981e-2 * (1LL << 32) + 0.5);
    const uint32_t c2 = (uint32_t)(4.9999838648363948e-1 * (1LL << 32) + 0.5);
    const uint32_t c3 = (uint32_t)(9.9999997159466147e-1 * (1LL << 32) + 0.5);
    // auxilliary constants
    const int32_t hpi_p23 = (int32_t)(3.141592653590 / 2 * (1LL << 23) + 0.5);
    const int32_t qpi_p23 = (int32_t)(3.141592653590 / 4 * (1LL << 23) + 0.5);
    const int32_t one_p23 = (int32_t)(1.0000000000000e+0 * (1LL << 23) + 0.5);
    uint32_t a, s, q, h, l, t, sn, cs;

    /* reduce range from [0, PI/2] to [0, PI/4] */
    t = (x > qpi_p23) ? (hpi_p23 - x) : x; // S8.23

    /* scale up argument for maximum precision in intermediate computation */
    a = t << 9; // U0.32

    /* pre-compute a**2 and a**4 */
    s = umul32hi (a, a); // U0.32
    q = umul32hi (s, s); // U0.32

    /* approximate sine on [0, PI/4] */
    h = s3 - umul32hi (s2, s); // U0.32
    l = umul32hi (s1 - umul32hi (s0, s), q); // U0.32
    sn = umul32hi (h + l, a); // U0.32

    /* approximate cosine on [0, PI/4] */
    h = c3 - umul32hi (c2, s); // U0.32
    l = umul32hi (c1 - umul32hi (c0, s), q); // U0.32
    cs = h + l; // U0.32

    /* round results to target precision */
    sn = ((sn + 256) >> 9); // S8.23
    cs = ((cs + 256) >> 9); // S8.23

    /* cosine result overflows U0.32 format for small arguments */
    cs = (t < 0xb50) ? one_p23 : cs; // S8.23

    /* map sine/cosine approximations based on quadrant */
    *sine   = (t != x) ? cs : sn; // S8.23
    *cosine = (t = x) ? sn : cs; // S8.23
}   

/*
  compute sine and cosine of argument in [0, PI/2]
  input and output in S8.23 format
  max err sine = 1.13173883e-7  max err cosine = 1.13158773e-7
  rms err sine = 4.30955921e-8  rms err cosine = 4.31472191e-8
  sin correctly rounded: 10844170 (82.30%)  
  cos correctly rounded: 10855609 (82.38%)

  Based on an approach by OllyW (http://www.olliw.eu/2014/fast-functions/, 
  retrieved 10/23/2020). We transform a = 2/PI*x-1/2, then we approximate
  sin ((2*PI*a + PI)/4 and cos ((2*PI*a + PI)/4. Except for sign flipping
  in the odd-power terms of the expansions the two series expansions match:

https://www.wolframalpha.com/input/?i=series++sin+%28%282*pi*a+%2B+pi%29%2F4%29
https://www.wolframalpha.com/input/?i=series++cos+%28%282*pi*a+%2B+pi%29%2F4%29

  This means we can sum the odd-power and the even-power terms seperately,
  then compute the sum and difference of those sums giving sine and cosine.
*/
void sincos_fixed_ollyw (int32_t x, int32_t *sine, int32_t *cosine)
{
    // minimax polynomial approximation for sin ((2*PI*a + PI)/4 on [-0.5, 0.5]
    const uint32_t c0 = (uint32_t)(7.0710676768794656e-1 * (1LL << 32) + 0.5);
    const uint32_t c1 = (uint32_t)((1.110721191857 -.25) * (1LL << 32) + 0.5);
    const uint32_t c2 = (uint32_t)(8.7235601339489222e-1 * (1LL << 32) + 0.5);
    const uint32_t c3 = (uint32_t)(4.5677902549505234e-1 * (1LL << 32) + 0.5);
    const uint32_t c4 = (uint32_t)(1.7932640877552330e-1 * (1LL << 32) + 0.5);
    const uint32_t c5 = (uint32_t)(5.6449491763487458e-2 * (1LL << 32) + 0.5);
    const uint32_t c6 = (uint32_t)(1.4444266213104129e-2 * (1LL << 32) + 0.5);
    const uint32_t c7 = (uint32_t)(3.4931597765535116e-3 * (1LL << 32) + 0.5);
    // auxiliary constants
    const uint32_t twoopi = (uint32_t)(2/3.1415926535898 * (1LL << 32) + 0.5);
    const uint32_t half_p31 = (uint32_t)(0.5000000000000 * (1LL << 31) + 0.5);
    const uint32_t quarter_p30 = (uint32_t)(0.2500000000 * (1LL << 30) + 0.5);
    uint32_t s, t, q, h, l;
    int32_t a, o, e, sn, cs;

    /* scale up argument for maximum precision in intermediate computation */
    t = (uint32_t)x << 8; // U1.31

    /* a = 2*PI*x - 0.5 */
    a = umul32hi (twoopi, t) - half_p31; // S0.31

    /* precompute a**2 and a**4 */
    s = (uint32_t)mul32hi (a, a) << 2; // U0.32
    q = umul32hi (s, s); // U0.32

    /* sum odd power terms; add in second portion of c1 (= 0.25) at the end */
    h = c1 - umul32hi (c3, s); // U0.32
    l = umul32hi ((c5 - umul32hi (c7, s)), q); // U0.32
    o = ((h + l) >> 2) + quarter_p30; // S1.30
    o = mul32hi (o, a); // S2.29

    /* sum even power terms */
    h = c0 - umul32hi (c2, s); // U0.32
    l = umul32hi ((c4 - umul32hi (c6, s)), q); // U0.32
    e = (h + l) >> 3; // S2.29 

    /* compute sine and cosine as sum and difference of odd / even terms */
    sn = e + o; // S2.29 sum -> sine 
    cs = e - o; // S2.29 difference -> cosine

    /* round results to target precision */
    sn = (sn + 32) >> 6; // S8.23
    cs = (cs + 32) >> 6; // S8.23

    *sine = sn;
    *cosine = cs;
}

double s8p23_to_double (int32_t a)
{
    return (double)a / (1LL << 23);
}

int32_t double_to_s8p23 (double a)
{
    return (int32_t)(a * (1LL << 23) + 0.5);
}

/* exhaustive test of S8.23 fixed-point sincos on [0,PI/2] */
int main (void)
{
    double errc, errs, maxerrs, maxerrc, errsqs, errsqc;
    int32_t arg, sin_correctly_rounded, cos_correctly_rounded;

#if VARIANT == SINCOS_OLLYW
    printf ("S8.23 fixed-point sincos OllyW variant\n");
#elif VARIANT == SINCOS_NJ
    printf ("S8.23 fixed-point sincos NJ variant\n");
#else // VARIANT
#error unsupported VARIANT
#endif // VARIANT

    maxerrs = 0; 
    maxerrc = 0;
    errsqs = 0;
    errsqc = 0;
    sin_correctly_rounded = 0;
    cos_correctly_rounded = 0;

    for (arg = 0; arg <= double_to_s8p23 (3.14159265358979 / 2); arg++) {
        double argf, refs, refc;
        int32_t sine, cosine, refsi, refci;
#if VARIANT == SINCOS_OLLYW
        sincos_fixed_ollyw (arg, &sine, &cosine);
#elif VARIANT == SINCOS_NJ
        sincos_fixed_nj (arg, &sine, &cosine);
#endif // VARIANT
        argf = s8p23_to_double (arg);
        refs = sin (argf);
        refc = cos (argf);
        refsi = double_to_s8p23 (refs);
        refci = double_to_s8p23 (refc);
        /* print function values near endpoints of interval */
        if ((arg < 5) || (arg > 0xc90fd5)) {
            printf ("arg=%08x  sin=%08x  cos=%08x\n", arg, sine, cosine);
        }
        if (sine == refsi) sin_correctly_rounded++;
        if (cosine == refci) cos_correctly_rounded++;
        errs = fabs (s8p23_to_double (sine) - refs);
        errc = fabs (s8p23_to_double (cosine) - refc);
        errsqs += errs * errs;
        errsqc += errc * errc;
        if (errs > maxerrs) maxerrs = errs;
        if (errc > maxerrc) maxerrc = errc;
    }
    printf ("max err sine = %15.8e  max err cosine = %15.8e\n", 
            maxerrs, maxerrc);
    printf ("rms err sine = %15.8e  rms err cosine = %15.8e\n", 
            sqrt (errsqs / arg), sqrt (errsqc / arg));
    printf ("sin correctly rounded: %d (%.2f%%)  cos correctly rounded: %d (%.2f%%)\n", 
            sin_correctly_rounded, 1.0 * sin_correctly_rounded / arg * 100,
            cos_correctly_rounded, 1.0 * cos_correctly_rounded / arg * 100);
    return EXIT_SUCCESS;
}

La salida del marco de pruebas adjunto debería ser esencialmente así:

S8.23 fixed-point sincos NJ variant
arg=00000000  sin=00000000  cos=00800000
arg=00000001  sin=00000001  cos=00800000
arg=00000002  sin=00000002  cos=00800000
arg=00000003  sin=00000003  cos=00800000
arg=00000004  sin=00000004  cos=00800000
arg=00c90fd6  sin=00800000  cos=00000005
arg=00c90fd7  sin=00800000  cos=00000004
arg=00c90fd8  sin=00800000  cos=00000003
arg=00c90fd9  sin=00800000  cos=00000002
arg=00c90fda  sin=00800000  cos=00000001
arg=00c90fdb  sin=00800000  cos=00000000
max err sine = 9.86237533e-008  max err cosine = 1.02729891e-007
rms err sine = 4.11141973e-008  rms err cosine = 4.11752018e-008
sin correctly rounded: 10961278 (83.19%)  cos correctly rounded: 11070113 (84.01%)

fixed-point sincos OllyW variant
arg=00000000  sin=00000000  cos=00800000
arg=00000001  sin=00000001  cos=00800000
arg=00000002  sin=00000002  cos=00800000
arg=00000003  sin=00000003  cos=00800000
arg=00000004  sin=00000004  cos=00800000
arg=00c90fd6  sin=00800000  cos=00000005
arg=00c90fd7  sin=00800000  cos=00000004
arg=00c90fd8  sin=00800000  cos=00000003
arg=00c90fd9  sin=00800000  cos=00000002
arg=00c90fda  sin=00800000  cos=00000001
arg=00c90fdb  sin=00800000  cos=00000000
max err sine = 1.13173883e-007  max err cosine = 1.13158773e-007
rms err sine = 4.30955921e-008  rms err cosine = 4.31472191e-008
sin correctly rounded: 10844170 (82.30%)  cos correctly rounded: 10855609 (82.38%)