Solución:
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, aquí hay una función sinusoidal que es igualmente rápida:
double FastSin(double x)
{
return x;
}
Esta respuesta en realidad no apesta, cuando x está cerca de cero. Para x pequeña, sin (x) es aproximadamente igual ax, porque x es el primer término de la expansión de Taylor de sin (x).
¿Qué, todavía no es lo suficientemente preciso para ti? 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 ciencias de la computación.
Debes 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” serán incorrectas.
La mayoría de las personas que hacen esta pregunta no comprenden la importancia de la compensaciones entre rendimiento y 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. Vea todas esas respuestas aleatorias debajo de la mía, que publican un montón de código fuente aleatorio sin comentar, sin documentar claramente el algoritmo utilizado y su exacto error máximo en todo el rango de entrada? “El error es aproximadamente una especie de murmullo, supongo”. Eso es estrictamente la liga de Bush. Si no sabe cómo calcular el PRECISO error máximo, a LLENO precisión, en su función de aproximación, a través del COMPLETO rango de las entradas … ¡entonces no sabes cómo escribir una función de aproximación!
Nadie usa la serie de Taylor por sí sola para aproximarse a los trascendentales en el software. A excepción de ciertos casos muy específicos, las series de Taylor generalmente se acercan al objetivo lentamente a través de rangos de entrada comunes.
Los algoritmos que usaban sus abuelos para calcular los trascendentales de manera eficiente, se conocen colectivamente como 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 rendimiento por precisión, incluida la que vinculé.
Ha habido muchas mejoras incrementales en los algoritmos CORDIC originales a lo largo de los años. Por ejemplo, el año pasado algunos investigadores de Japón publicaron un artículo sobre un CORDIC mejorado con mejores ángulos de rotación, lo que reduce las operaciones necesarias.
Si tiene multiplicadores de hardware por ahí (y es casi seguro que los tenga), 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, en 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 en un rango enorme, ignorando la simetría natural de la función seno y simplemente resolviendo el problema de aproximación arrojándole más coeficientes. Y aquí hay un ejemplo de cómo estimar una función sinusoidal dentro de 5 ULP. ¿No sabes qué es una ULP? Debería.
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 procesando señales digitales, los polinomios de Chebyshev le brindan un efecto de difuminado económico y predecible “gratis”.
Si desea encontrar sus propios coeficientes polinomiales de Chebyshev en un rango específico, muchas bibliotecas matemáticas llaman al proceso de encontrar esos coeficientes “ajuste de 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. Por lo general, cuando alguien crea un algoritmo “increíblemente nuevo” para hacer raíces cuadradas, es simplemente Newton-Raphson disfrazado.
Los polinomios de Newton-Raphson, CORDIC y Chebyshev le permiten intercambiar velocidad por precisión, por lo que la respuesta puede ser tan imprecisa como desee.
Por último, cuando haya terminado con todas sus elegantes evaluaciones comparativas y microoptimización, asegúrese de que su versión “rápida” sea en realidad más rápida que la versión de la biblioteca. Aquí hay una implementación de biblioteca típica de fsin () limitada al dominio de -pi / 4 a pi / 4. Y no es tan condenadamente lento.
Una última advertencia para usted: es casi seguro que esté utilizando las matemáticas IEEE-754 para realizar sus estimaciones, y cada vez que realice las matemáticas IEEE-754 con un montón de multiplicaciones, algunas decisiones de ingeniería oscuras tomadas hace décadas volverán a acechar. usted, en forma de errores de redondeo. ¡Y esos errores comienzan pequeños, pero se hacen cada vez más grandes y más grandes y más GRANDES! En algún momento de su vida, lea “Lo que todo científico informático debería saber sobre los números de coma flotante” y tenga la cantidad adecuada de miedo. Tenga en cuenta que si comienza a escribir sus propias funciones trascendentales, deberá comparar y medir el error REAL debido al redondeo de punto flotante, no solo el error teórico máximo. Ésta no es una preocupación teórica; La configuración de compilación de “matemáticas rápidas” me ha mordido en el trasero, en más de un proyecto.
tl: dr; vaya a google “aproximación de seno” o “aproximación de coseno” o “aproximación de raíz cuadrada” o “teoría de aproximación”.
Primero, la serie Taylor NO es la mejor / más rápida forma de implementar seno / 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 velocidad de manera más eficiente. Además, este problema ya se ha discutido ampliamente en StackOverflow. He aquí solo un ejemplo.
En segundo lugar, la gran diferencia que ve entre PCS antiguo / nuevo se debe al hecho de que la arquitectura Intel moderna tiene un código de ensamblaje explícito para calcular funciones trigonométricas elementales. Es bastante difícil vencerlos en velocidad de ejecución.
Por último, hablemos del código de su antigua PC. Verifique la implementación de la biblioteca científica gsl gnu (o recetas numéricas), y verá que básicamente usan 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é detalles de implementación aquí porque ya hay muy buenas respuestas publicadas en StackOverflow. Mira este, por ejemplo. Simplemente modifique el número de términos de esta serie para cambiar el equilibrio entre precisión / velocidad.
Para este tipo de problema, si desea obtener detalles de implementación de alguna función especial o método numérico, debería echar un vistazo al 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 al incluir indicadores de optimización de punto flotante agresivos en gcc / icc. Esto disminuirá la precisión, pero parece que eso es exactamente lo que desea.
EDIT2: Puede intentar hacer una cuadrícula de sin gruesa 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 pérdida de caché)!
Para la raíz cuadrada, existe un enfoque llamado desplazamiento de bits.
Un número flotante definido por IEEE-754 utiliza algunos bits específicos para representar tiempos de descripción de múltiples basados en 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;
}
Eso es un tiempo constante para calcular la raíz cuadrada.