simbolo raiz raices radicacion programa enesima dev cuadrada c++ math optimization trigonometry

raiz - La implementación más rápida de seno, coseno y raíz cuadrada en C++(no necesita ser muy precisa)



raiz enesima en c++ (11)

He estado buscando en Google la pregunta durante la última hora, pero solo hay puntos para Taylor Series o algún código de ejemplo que es demasiado lento o no se compila en absoluto. Bueno, la mayoría de las respuestas que he encontrado en Google es "Google, ya está preguntado", pero lamentablemente no es ...

Estoy perfilando mi juego en el Pentium 4 de gama baja y descubrí que ~ 85% del tiempo de ejecución se desperdicia en el cálculo de sinus, cosinus y raíz cuadrada (de la biblioteca estándar de C ++ en Visual Studio), y esto parece ser muy dependiente de la CPU ( en mi I7 las mismas funciones obtuvieron solo el 5% del tiempo de ejecución, y el juego es waaaaaaaaaay más rápido). No puedo optimizar estas tres funciones, ni calcular el seno y el coseno en una sola pasada (hay interdependientes), pero no necesito resultados demasiado precisos para mi simulación, por lo que puedo vivir con una aproximación más rápida.

Entonces, la pregunta: ¿Cuál es la forma más rápida de calcular el seno, el coseno y la raíz cuadrada para la flotación en C ++?

EDITAR La tabla de búsqueda es más dolorosa, ya que la pérdida de caché resultante es mucho más costosa en la CPU moderna que la serie Taylor. Los CPU son tan rápidos en estos días, y el caché no lo es.

Cometí un error, pensé que necesitaba calcular varios factoriales para la serie de Taylor, y ahora veo que se pueden implementar como constantes.

Así que la pregunta actualizada: ¿hay alguna optimización rápida para la raíz cuadrada también?

EDIT2

Estoy usando la raíz cuadrada para calcular la distancia, no la normalización; no puedo usar el algoritmo de la raíz cuadrada inversa rápida (como se señaló en el comentario: http://en.wikipedia.org/wiki/Fast_inverse_square_root

EDITAR3

Tampoco puedo operar en distancias cuadradas, necesito una distancia exacta para los cálculos

Aquí está la función sinusoidal más rápida posible garantizada en C ++:

double FastSin(double x) { return 0; }

Oh, querías una precisión mejor que | 1.0 |? Bueno sigue leyendo.

Los ingenieros de la década de 1970 hicieron algunos descubrimientos fantásticos en este campo, pero los nuevos programadores simplemente desconocen que estos métodos existen, porque no se enseñan como parte de los planes de estudio estándar de informática.

Debe comenzar por comprender que no hay una implementación "perfecta" de estas funciones para todas las aplicaciones. Por lo tanto, se garantiza que las respuestas superficiales a preguntas como "cuál es el más rápido" están equivocadas.

La mayoría de las personas que hacen esta pregunta no comprenden la importancia de las compensaciones entre el rendimiento y la precisión . En particular, tendrá que tomar algunas decisiones con respecto a la precisión de los cálculos antes de hacer cualquier otra cosa. ¿Cuánto error puedes tolerar en el resultado? 10 ^ -4? 10 ^ -16?

A menos que pueda cuantificar el error en cualquier método, no lo use.

Nadie usa la serie de Taylor solo para aproximar los trascendentales en el software. A excepción de ciertos casos altamente específicos, las series de Taylor generalmente se acercan al objetivo lentamente a través de rangos de entrada comunes.

Los algoritmos que utilizaron sus abuelos para calcular los trascendentales de manera eficiente, se denominan colectivamente CORDIC y eran lo suficientemente simples como para ser implementados en hardware. Aquí hay una implementación CORDIC bien documentada en C. Las implementaciones CORDIC, por lo general, requieren una tabla de búsqueda muy pequeña, pero la mayoría de las implementaciones ni siquiera requieren que haya un multiplicador de hardware disponible. La mayoría de las implementaciones de CORDIC le permiten intercambiar el rendimiento por la precisión, incluido el que vinculé.

Se han producido muchas mejoras incrementales en los algoritmos CORDIC originales a lo largo de los años. Por ejemplo, el año pasado, algunos investigadores en Japón publicaron un article sobre un CORDIC mejorado con mejores ángulos de rotación, lo que reduce las operaciones requeridas.

Si tiene multiplicadores de hardware alrededor (y es casi seguro que sí), o si no puede pagar una tabla de búsqueda como la que requiere CORDIC, siempre puede usar un polinomio de Chebyshev para hacer lo mismo. Los polinomios de Chebyshev requieren multiplicaciones, pero esto rara vez es un problema en el hardware moderno. Nos gustan los polinomios de Chebyshev porque tienen errores máximos altamente predecibles para una aproximación dada . El máximo del último término en un polinomio de Chebyshev, a través de su rango de entrada, limita el error en el resultado. Y este error se reduce a medida que aumenta el número de términos. Aquí hay un ejemplo de un polinomio de Chebyshev que da una aproximación sinusoidal a través de un amplio rango, ignorando la simetría natural de la función seno y simplemente resolviendo el problema de aproximación arrojándole más coeficientes.

También nos gustan los polinomios de Chebyshev porque el error en la aproximación se distribuye por igual en el rango de salidas. Si está escribiendo complementos de audio o haciendo procesamiento de señal digital, los polinomios de Chebyshev le brindan un efecto de interpolación barato y predecible "de forma gratuita".

Si desea encontrar sus propios coeficientes polinómicos de Chebyshev en un rango específico, muchas bibliotecas de matemáticas consideran que el proceso de encontrar esos coeficientes "se ajusta a Chebyshev " o algo así.

Las raíces cuadradas, entonces como ahora, generalmente se calculan con alguna variante del algoritmo de Newton-Raphson , generalmente con un número fijo de iteraciones. Generalmente, cuando alguien crea un algoritmo "nuevo y sorprendente" para hacer raíces cuadradas, es simplemente Newton-Raphson disfrazado.

Los polinomios de Newton-Raphson, CORDIC y Chebyshev te permiten cambiar la velocidad por la precisión, por lo que la respuesta puede ser tan imprecisa como quieras.

Por último, cuando haya finalizado todas sus sofisticadas evaluaciones comparativas y micro optimización, asegúrese de que su versión "rápida" sea realmente más rápida que la versión de la biblioteca. Aquí hay una implementación de biblioteca típica de fsin () limitada en el dominio desde -pi / 4 a pi / 4. Y no es tan condenadamente lento.

Hay personas que han dedicado sus vidas a resolver estos problemas de manera eficiente y han producido algunos resultados fascinantes. Cuando estés listo para unirte a la vieja escuela, recoge una copia de Recetas numéricas .

tl: dr; Vaya a google "aproximación senoidal" o "aproximación coseno" o "aproximación de raíz cuadrada" o " teoría de aproximación ".

Aquí tienes gente, come y dame todos tus créditos, por favor y gracias:

void sincos_fast(float x, float *pS, float *pC){ float cosOff4LUT[] = { 0x1.000000p+00, 0x1.6A09E6p-01, 0x0.000000p+00, -0x1.6A09E6p-01, -0x1.000000p+00, -0x1.6A09E6p-01, 0x0.000000p+00, 0x1.6A09E6p-01 }; int m, ms, mc; float xI, xR, xR2; float c, s, cy, sy; // Cody & Waite''s range reduction Algorithm, [-pi/4, pi/4] xI = floorf(x * 0x1.45F306p+00 + 0.5); xR = (x - xI * 0x1.920000p-01) - xI*0x1.FB5444p-13; m = (int) xI; xR2 = xR*xR; // Find cosine & sine index for angle offsets indices mc = ( m ) & 0x7; // two''s complement permits upper modulus for negative numbers =P ms = (m + 6) & 0x7; // two''s complement permits upper modulus for negative numbers =P, note phase correction for sine. // Find cosine & sine cy = cosOff4LUT[mc]; // Load angle offset neighborhood cosine value sy = cosOff4LUT[ms]; // Load angle offset neighborhood sine value c = 0xf.ff79fp-4 + x2 * (-0x7.e58e9p-4); // TOL = 1.2786e-4 // c = 0xf.ffffdp-4 + xR2 * (-0x7.ffebep-4 + xR2 * 0xa.956a9p-8); // TOL = 1.7882e-7 s = xR * (0xf.ffbf7p-4 + x2 * (-0x2.a41d0cp-4)); // TOL = 4.835251e-6 // s = xR * (0xf.fffffp-4 + xR2 * (-0x2.aaa65cp-4 + xR2 * 0x2.1ea25p-8)); // TOL = 1.1841e-8 *pC = c*cy - s*sy; *pS = c*sy + s*cy; } float sqrt_fast(float x){ union {float f; int i; } X, Y; float ScOff; uint8_t e; Xf = x; e = (Xi >> 23); // f.SFPbits.e; if(x <= 0) return(0.0f); ScOff = ((e & 1) != 0) ? 1.0f : 0x1.6a09e6p0; // NOTE: If exp=EVEN, b/c (exp-127) a (EVEN - ODD) := ODD; but a (ODD - ODD) := EVEN!! e = ((e + 127) >> 1); // NOTE: If exp=ODD, b/c (exp-127) then flr((exp-127)/2) Xi = (Xi & ((1uL << 23) - 1)) | (0x7F << 23); // Mask mantissa, force exponent to zero. Yi = (((uint32_t) e) << 23); // Error grows with square root of the exponent. Unfortunately no work around like inverse square root... :( // Yf *= ScOff * (0x9.5f61ap-4 + Xf*(0x6.a09e68p-4)); // Error = +-1.78e-2 * 2^(flr(log2(x)/2)) // Yf *= ScOff * (0x7.2181d8p-4 + Xf*(0xa.05406p-4 + Xf*(-0x1.23a14cp-4))); // Error = +-7.64e-5 * 2^(flr(log2(x)/2)) // Yf *= ScOff * (0x5.f10e7p-4 + Xf*(0xc.8f2p-4 +Xf*(-0x2.e41a4cp-4 + Xf*(0x6.441e6p-8)))); // Error = 8.21e-5 * 2^(flr(log2(x)/2)) // Yf *= ScOff * (0x5.32eb88p-4 + Xf*(0xe.abbf5p-4 + Xf*(-0x5.18ee2p-4 + Xf*(0x1.655efp-4 + Xf*(-0x2.b11518p-8))))); // Error = +-9.92e-6 * 2^(flr(log2(x)/2)) // Yf *= ScOff * (0x4.adde5p-4 + Xf*(0x1.08448cp0 + Xf*(-0x7.ae1248p-4 + Xf*(0x3.2cf7a8p-4 + Xf*(-0xc.5c1e2p-8 + Xf*(0x1.4b6dp-8)))))); // Error = +-1.38e-6 * 2^(flr(log2(x)/2)) // Yf *= ScOff * (0x4.4a17fp-4 + Xf*(0x1.22d44p0 + Xf*(-0xa.972e8p-4 + Xf*(0x5.dd53fp-4 + Xf*(-0x2.273c08p-4 + Xf*(0x7.466cb8p-8 + Xf*(-0xa.ac00ep-12))))))); // Error = +-2.9e-7 * 2^(flr(log2(x)/2)) Yf *= ScOff * (0x3.fbb3e8p-4 + Xf*(0x1.3b2a3cp0 + Xf*(-0xd.cbb39p-4 + Xf*(0x9.9444ep-4 + Xf*(-0x4.b5ea38p-4 + Xf*(0x1.802f9ep-4 + Xf*(-0x4.6f0adp-8 + Xf*(0x5.c24a28p-12 )))))))); // Error = +-2.7e-6 * 2^(flr(log2(x)/2)) return(Yf); }

Basada en la idea de http://forum.devmaster.net/t/fast-and-accurate-sine-cosine/9648 y alguna reescritura manual para mejorar el rendimiento en un micro benchmark, terminé con la siguiente implementación de coseno que es se utiliza en una simulación física de HPC que se ve obstaculizada por repetidas llamadas a cos en un gran número de espacio. Es lo suficientemente preciso y mucho más rápido que una tabla de búsqueda, sobre todo no se requiere división.

template<typename T> inline T cos(T x) noexcept { constexpr T tp = 1./(2.*M_PI); x *= tp; x -= T(.25) + std::floor(x + T(.25)); x *= T(16.) * (std::abs(x) - T(.5)); #if EXTRA_PRECISION x += T(.225) * x * (std::abs(x) - T(1.)); #endif return x; }

El compilador de Intel, al menos, también es lo suficientemente inteligente como para vectorizar esta función cuando se usa en un bucle.

Si se define EXTRA_PRECISION, el error máximo es aproximadamente 0.00109 para el rango -π a π, asumiendo que T es el double como se define generalmente en la mayoría de las implementaciones de C ++. De lo contrario, el error máximo es aproximadamente 0.056 para el mismo rango.

Esta es una implementación del método de la serie Taylor que se dio anteriormente en la respuesta de akellehe .

unsigned int Math::SIN_LOOP = 15; unsigned int Math::COS_LOOP = 15; // sin(x) = x - x^3/3! + x^5/5! - x^7/7! + ... template <class T> T Math::sin(T x) { T Sum = 0; T Power = x; T Sign = 1; const T x2 = x * x; T Fact = 1.0; for (unsigned int i=1; i<SIN_LOOP; i+=2) { Sum += Sign * Power / Fact; Power *= x2; Fact *= (i + 1) * (i + 2); Sign *= -1.0; } return Sum; } // cos(x) = 1 - x^2/2! + x^4/4! - x^6/6! + ... template <class T> T Math::cos(T x) { T Sum = x; T Power = x; T Sign = 1.0; const T x2 = x * x; T Fact = 1.0; for (unsigned int i=3; i<COS_LOOP; i+=2) { Power *= x2; Fact *= i * (i - 1); Sign *= -1.0; Sum += Sign * Power / Fact; } return Sum; }

La forma más rápida es pre-calcular los valores y usar una tabla como en este ejemplo:

Crear tabla de búsqueda de seno en C ++

PERO si insiste en computar en tiempo de ejecución, puede usar la expansión de seno o coseno de la serie Taylor ...

Para más información sobre la serie de Taylor ... http://en.wikipedia.org/wiki/Taylor_series

Una de las claves para que esto funcione bien es pre-computar los factoriales y truncar en un número razonable de términos. Los factoriales crecen en el denominador muy rápidamente, por lo que no necesita tener más de unos pocos términos.

Además ... No multipliques tu x ^ n desde el principio cada vez ... por ejemplo, multiplica x ^ 3 por x otras dos veces, luego por otras dos para calcular los exponentes.

Más de 100000000 de prueba, milianw answer es 2 veces más lenta que la implementación de std :: cos. Sin embargo, puedes administrar para ejecutarlo más rápido siguiendo los siguientes pasos:

-> usar flotador

-> no use floor pero static_cast

-> No usar abdominales sino condicionales ternarios.

-> use #define constante para la división

-> usar macro para evitar la llamada de función

// 1 / (2 * PI) #define FPII 0.159154943091895 //PI / 2 #define PI2 1.570796326794896619 #define _cos(x) x *= FPII;/ x -= .25f + static_cast<int>(x + .25f) - 1;/ x *= 16.f * ((x >= 0 ? x : -x) - .5f); #define _sin(x) x -= PI2; _cos(x);

Más de 100000000 llamadas a std :: cos y _cos (x), std :: cos se ejecutan en ~ 14s vs ~ 3s para _cos (x) (un poco más para _sin (x))

Para la raíz cuadrada, hay un enfoque llamado cambio de bits.

Un número flotante definido por IEEE-754 está utilizando algunos bits determinados que representan tiempos de múltiples bases 2. Algunos bits son para representar el valor base.

float squareRoot(float x) { unsigned int i = *(unsigned int*) &x; // adjust bias i += 127 << 23; // approximation of square root i >>= 1; return *(float*) &i; }

Es un tiempo constante calculando la raíz cuadrada.

Para x86, las instrucciones de raíz cuadrada de hardware FP son rápidas ( sqrtss es sqrt Scalar Single-precision). La precisión simple es más rápida que la precisión doble, por lo que definitivamente use float lugar de double para el código donde puede permitirse usar menos precisión.

Para el código de 32 bits, por lo general necesita opciones de compilador para hacer FP FP con instrucciones SSE, en lugar de x87. (por ejemplo, -mfpmath=sse )

Para que las funciones sqrt() o sqrtf() C sqrtf() en línea como sqrtsd o sqrtss , debe compilar con -fno-math-errno . Tener funciones matemáticas configuradas errno en NaN generalmente se considera un error de diseño, pero el estándar lo requiere. Sin esa opción, gcc lo alinea, pero luego hace una comparación + rama para ver si el resultado fue NaN, y si es así, llama a la función de biblioteca para que pueda configurar errno . Si su programa no comprueba errno después de las funciones matemáticas, no hay peligro en usar -fno-math-errno .

No necesita ninguna de las partes "inseguras" de -ffast-math para obtener sqrt y algunas otras funciones para integrarse mejor o en absoluto, pero -ffast-math puede hacer una gran diferencia (por ejemplo, permitir que el compilador se auto-vectorice) en los casos en que eso cambia el resultado, porque la matemática FP no es asociativa .

por ejemplo, con gcc6.3 compilando float foo(float a){ return sqrtf(a); } float foo(float a){ return sqrtf(a); }

foo: # with -O3 -fno-math-errno. sqrtss xmm0, xmm0 ret

foo: # with just -O3 pxor xmm2, xmm2 # clang just checks for NaN, instead of comparing against zero. sqrtss xmm1, xmm0 ucomiss xmm2, xmm0 ja .L8 # take the slow path if 0.0 > a movaps xmm0, xmm1 ret .L8: # errno-setting path sub rsp, 24 movss DWORD PTR [rsp+12], xmm1 # store the sqrtss result because the x86-64 SysV ABI has no call-preserved xmm regs. call sqrtf # call sqrtf just to set errno movss xmm1, DWORD PTR [rsp+12] add rsp, 24 movaps xmm0, xmm1 # extra mov because gcc reloaded into the wrong register. ret

El código de gcc para el caso de NaN parece demasiado complicado; ¡Ni siquiera usa el valor de retorno de sqrtf ! De todos modos, este es el tipo de desorden que realmente obtienes sin -fno-math-errno , para cada sitio de llamadas sqrtf() en tu programa. En su mayoría, solo se trata de código inflado, y ninguno de los bloques .L8 se ejecutarán cuando se tome el cuadrado de un número> = 0.0, pero aún hay varias instrucciones adicionales en la ruta rápida.

Si sabe que su entrada a sqrt no es cero, puede usar la instrucción recíproca de sqrt rápida, pero muy aproximada, rsqrtps (o rsqrtss para la versión escalar). Una iteración de Newton-Raphson lo lleva a casi la misma precisión que la instrucción sqrt precisión simple del hardware, pero no del todo.

sqrt(x) = x * 1/sqrt(x) , para x!=0 , por lo que puede calcular un sqrt con rsqrt y una multiplicación. Ambos son rápidos, incluso en P4 (¿era eso aún relevante en 2013)?

En P4, puede valer la pena usar la iteración rsqrt + Newton para reemplazar un solo sqrt , incluso si no es necesario dividir nada por él.

Vea también una respuesta que escribí recientemente sobre el manejo de ceros al calcular sqrt(x) como x*rsqrt(x) , con una iteración de Newton. Incluí una discusión sobre el error de redondeo si desea convertir el valor de PF en un entero, y enlaces a otras preguntas relevantes.

P4:

  • sqrtss : 23c de latencia, no canalizado
  • sqrtsd : 38c de latencia, no canalizado
  • fsqrt (x87): 43c de latencia, no canalizado

  • rsqrtss / mulss : 4c + 6c de latencia. Posiblemente uno por 3c de rendimiento, ya que aparentemente no necesitan la misma unidad de ejecución (mmx vs. fp).

  • Las versiones empaquetadas SIMD son algo más lentas

Skylake:

  • sqrtss / sqrtps : 12c de latencia, una por rendimiento de 3c
  • sqrtsd / sqrtpd : latencia de 15-16c, una por rendimiento de 4-6c
  • fsqrt (x87): 14-21cc de latencia, una por 4-7c de rendimiento
  • rsqrtss / mulss : 4c + 4c de latencia. Uno por 1c de rendimiento.
  • Las versiones vectoriales de SIMD 128b tienen la misma velocidad. Las versiones vectoriales de 256b tienen una latencia un poco mayor, casi la mitad del rendimiento. La versión rsqrtss tiene un rendimiento completo para vectores de 256b.

Con una iteración de Newton, la versión rsqrt no es mucho más rápida.

Números de las pruebas experimentales de Agner Fog . Consulte sus guías de microarca para comprender qué hace que el código se ejecute rápido o lento. También vea los enlaces en la etiqueta wiki x86 .

IDK como mejor calcular el pecado / cos. He leído que el hardware fsin / fcos (y los únicos fsincos un poco más fsincos que hacen las dos cosas a la vez) no son la forma más rápida, pero lo que es IDK.

Primero, la serie de Taylor NO es la manera mejor / más rápida de implementar sine / cos. Tampoco es la forma en que las bibliotecas profesionales implementan estas funciones trigonométricas, y conocer la mejor implementación numérica le permite ajustar la precisión para obtener una velocidad más eficiente. Además, este problema ya se ha tratado ampliamente en . Aquí es sólo un ejemplo .

En segundo lugar, la gran diferencia que se ve entre PCS antiguos / nuevos se debe al hecho de que la arquitectura moderna de Intel tiene un código de ensamblaje explícito para calcular funciones trigonométricas elementales. Es bastante difícil vencerlos en velocidad de ejecución.

Finalmente, hablemos del código en tu vieja PC. Verifique la implementación de la biblioteca científica (o recetas numéricas) de gsl gnu , y verá que básicamente utilizan la fórmula de aproximación de Chebyshev.

La aproximación de Chebyshev converge más rápido, por lo que necesita evaluar menos términos. No escribiré los detalles de implementación aquí porque ya hay respuestas muy buenas publicadas en . Marque este por ejemplo . Solo modifique la cantidad de términos en esta serie para cambiar el equilibrio entre precisión / velocidad.

Por cierto: la regla 0 para este tipo de problema: si desea detalles de implementación de alguna función especial o método numérico, debe consultar el código GSL antes de realizar cualquier otra acción. GSL es la biblioteca numérica ESTÁNDAR.

EDITAR: puede mejorar el tiempo de ejecución incluyendo indicadores agresivos de optimización de punto flotante en gcc / icc. Esto disminuirá la precisión, pero parece que eso es exactamente lo que quieres.

EDIT2: puede intentar hacer una cuadrícula de pecado grueso y usar la rutina gsl (gsl_interp_cspline_periodic para spline con condiciones periódicas) para spline esa tabla (la spline reducirá los errores en comparación con una interpolación lineal => necesita menos puntos en su tabla = > menos falta de caché)!

QT tiene implementaciones rápidas de seno (qFastSin) y coseno (qFastCos) que utiliza la tabla de consulta con interpolación. Lo estoy usando en mi código y son más rápidos que std: sin / cos y lo suficientemente precisos para lo que necesito (error ~ = 0.01%, supongo):

https://code.woboq.org/qt5/qtbase/src/corelib/kernel/qmath.h.html#_Z8qFastSind

void sincos_fast(float x, float *pS, float *pC){ float cosOff4LUT[] = { 0x1.000000p+00, 0x1.6A09E6p-01, 0x0.000000p+00, -0x1.6A09E6p-01, -0x1.000000p+00, -0x1.6A09E6p-01, 0x0.000000p+00, 0x1.6A09E6p-01 }; int m, ms, mc; float xI, xR, xR2; float c, s, cy, sy; // Cody & Waite''s range reduction Algorithm, [-pi/4, pi/4] xI = floorf(x * 0x1.45F306p+00 + 0.5); xR = (x - xI * 0x1.920000p-01) - xI*0x1.FB5444p-13; m = (int) xI; xR2 = xR*xR; // Find cosine & sine index for angle offsets indices mc = ( m ) & 0x7; // two''s complement permits upper modulus for negative numbers =P ms = (m + 6) & 0x7; // two''s complement permits upper modulus for negative numbers =P, note phase correction for sine. // Find cosine & sine cy = cosOff4LUT[mc]; // Load angle offset neighborhood cosine value sy = cosOff4LUT[ms]; // Load angle offset neighborhood sine value c = 0xf.ff79fp-4 + x2 * (-0x7.e58e9p-4); // TOL = 1.2786e-4 // c = 0xf.ffffdp-4 + xR2 * (-0x7.ffebep-4 + xR2 * 0xa.956a9p-8); // TOL = 1.7882e-7 s = xR * (0xf.ffbf7p-4 + x2 * (-0x2.a41d0cp-4)); // TOL = 4.835251e-6 // s = xR * (0xf.fffffp-4 + xR2 * (-0x2.aaa65cp-4 + xR2 * 0x2.1ea25p-8)); // TOL = 1.1841e-8 *pC = c*cy - s*sy; *pS = c*sy + s*cy; } float sqrt_fast(float x){ union {float f; int i; } X, Y; float ScOff; uint8_t e; X.f = x; e = (X.i >> 23); // f.SFPbits.e; if(x <= 0) return(0.0f); ScOff = ((e & 1) != 0) ? 1.0f : 0x1.6a09e6p0; // NOTE: If exp=EVEN, b/c (exp-127) a (EVEN - ODD) := ODD; but a (ODD - ODD) := EVEN!! e = ((e + 127) >> 1); // NOTE: If exp=ODD, b/c (exp-127) then flr((exp-127)/2) X.i = (X.i & ((1uL << 23) - 1)) | (0x7F << 23); // Mask mantissa, force exponent to zero. Y.i = (((uint32_t) e) << 23); // Error grows with square root of the exponent. Unfortunately no work around like inverse square root... :( // Y.f *= ScOff * (0x9.5f61ap-4 + X.f*(0x6.a09e68p-4)); // Error = +-1.78e-2 * 2^(flr(log2(x)/2)) // Y.f *= ScOff * (0x7.2181d8p-4 + X.f*(0xa.05406p-4 + X.f*(-0x1.23a14cp-4))); // Error = +-7.64e-5 * 2^(flr(log2(x)/2)) // Y.f *= ScOff * (0x5.f10e7p-4 + X.f*(0xc.8f2p-4 +X.f*(-0x2.e41a4cp-4 + X.f*(0x6.441e6p-8)))); // Error = 8.21e-5 * 2^(flr(log2(x)/2)) // Y.f *= ScOff * (0x5.32eb88p-4 + X.f*(0xe.abbf5p-4 + X.f*(-0x5.18ee2p-4 + X.f*(0x1.655efp-4 + X.f*(-0x2.b11518p-8))))); // Error = +-9.92e-6 * 2^(flr(log2(x)/2)) // Y.f *= ScOff * (0x4.adde5p-4 + X.f*(0x1.08448cp0 + X.f*(-0x7.ae1248p-4 + X.f*(0x3.2cf7a8p-4 + X.f*(-0xc.5c1e2p-8 + X.f*(0x1.4b6dp-8)))))); // Error = +-1.38e-6 * 2^(flr(log2(x)/2)) // Y.f *= ScOff * (0x4.4a17fp-4 + X.f*(0x1.22d44p0 + X.f*(-0xa.972e8p-4 + X.f*(0x5.dd53fp-4 + X.f*(-0x2.273c08p-4 + X.f*(0x7.466cb8p-8 + X.f*(-0xa.ac00ep-12))))))); // Error = +-2.9e-7 * 2^(flr(log2(x)/2)) Y.f *= ScOff * (0x3.fbb3e8p-4 + X.f*(0x1.3b2a3cp0 + X.f*(-0xd.cbb39p-4 + X.f*(0x9.9444ep-4 + X.f*(-0x4.b5ea38p-4 + X.f*(0x1.802f9ep-4 + X.f*(-0x4.6f0adp-8 + X.f*(0x5.c24a28p-12 )))))))); // Error = +-2.7e-6 * 2^(flr(log2(x)/2)) return(Y.f); }

La LUT y la licencia se pueden encontrar aquí: https://code.woboq.org/qt5/qtbase/src/corelib/kernel/qmath.cpp.html#qt_sine_table

Utilizo el siguiente código para calcular funciones trigonométricas con una precisión cuádruple. La constante N determina el número de bits de precisión requeridos (por ejemplo, N = 26 dará una precisión de precisión única). Dependiendo de la precisión deseada, el almacenamiento precomputado puede ser pequeño y cabrá en el caché. Solo requiere operaciones de suma y multiplicación y también es muy fácil de vectorizar.

El algoritmo calcula previamente los valores de sin y cos para 0.5 ^ i, i = 1, ..., N. Luego, podemos combinar estos valores precomputados, para calcular sin y cos para cualquier ángulo hasta una resolución de 0.5 ^ N

template <class QuadReal_t> QuadReal_t sin(const QuadReal_t a){ const int N=128; static std::vector<QuadReal_t> theta; static std::vector<QuadReal_t> sinval; static std::vector<QuadReal_t> cosval; if(theta.size()==0){ #pragma omp critical (QUAD_SIN) if(theta.size()==0){ theta.resize(N); sinval.resize(N); cosval.resize(N); QuadReal_t t=1.0; for(int i=0;i<N;i++){ theta[i]=t; t=t*0.5; } sinval[N-1]=theta[N-1]; cosval[N-1]=1.0-sinval[N-1]*sinval[N-1]/2; for(int i=N-2;i>=0;i--){ sinval[i]=2.0*sinval[i+1]*cosval[i+1]; cosval[i]=sqrt(1.0-sinval[i]*sinval[i]); } } } QuadReal_t t=(a<0.0?-a:a); QuadReal_t sval=0.0; QuadReal_t cval=1.0; for(int i=0;i<N;i++){ while(theta[i]<=t){ QuadReal_t sval_=sval*cosval[i]+cval*sinval[i]; QuadReal_t cval_=cval*cosval[i]-sval*sinval[i]; sval=sval_; cval=cval_; t=t-theta[i]; } } return (a<0.0?-sval:sval); }