He estado estudiando detenidamente los desensamblajes de .NET y el código fuente de GCC, pero parece que no puedo encontrar en ninguna parte la implementación real sin()
y otras funciones matemáticas ... siempre parecen estar haciendo referencia a otra cosa.
¿Alguien puede ayudarme a encontrarlos? Siento que es poco probable que TODO el hardware en el que se ejecutará C admita funciones trigonométricas en el hardware, por lo que debe haber un algoritmo de software en algún lugar , ¿verdad?
Soy consciente de varias maneras en que se pueden calcular las funciones , y he escrito mis propias rutinas para calcular funciones usando la serie taylor por diversión. Tengo curiosidad por saber cuán reales son los lenguajes de producción, ya que todas mis implementaciones son siempre varios órdenes de magnitud más lentas, aunque creo que mis algoritmos son bastante inteligentes (obviamente no lo son).
fuente
Respuestas:
En GNU libm, la implementación de
sin
depende del sistema. Por lo tanto, puede encontrar la implementación, para cada plataforma, en algún lugar del subdirectorio apropiado de sysdeps .Un directorio incluye una implementación en C, aportada por IBM. Desde octubre de 2011, este es el código que realmente se ejecuta cuando se llama
sin()
a un sistema Linux x86-64 típico. Aparentemente es más rápido que lasfsin
instrucciones de montaje. Código fuente: sysdeps / ieee754 / dbl-64 / s_sin.c , busque__sin (double x)
.Este código es muy complejo. Ningún algoritmo de software es tan rápido como sea posible y también preciso en todo el rango de valores de x , por lo que la biblioteca implementa varios algoritmos diferentes, y su primer trabajo es mirar xy decidir qué algoritmo usar.
Cuando x está muy muy cerca de 0,
sin(x) == x
es la respuesta correcta.Un poco más lejos,
sin(x)
utiliza la serie familiar de Taylor. Sin embargo, esto solo es exacto cerca de 0, así que ...Cuando el ángulo es mayor de aproximadamente 7 °, se usa un algoritmo diferente, que calcula las aproximaciones de la serie Taylor para sin (x) y cos (x), y luego usa valores de una tabla calculada previamente para refinar la aproximación.
Cuando | x | > 2, ninguno de los algoritmos anteriores funcionaría, por lo que el código comienza calculando algún valor más cercano a 0 que puede ser alimentado
sin
ocos
no.Hay otra rama para tratar con x siendo un NaN o infinito.
Este código usa algunos hacks numéricos que nunca antes había visto, aunque por lo que sé, podrían ser conocidos entre los expertos en coma flotante. A veces, algunas líneas de código tomarían varios párrafos para explicar. Por ejemplo, estas dos líneas.
se usan (a veces) para reducir x a un valor cercano a 0 que difiere de x en un múltiplo de π / 2, específicamente
xn
× π / 2. La forma en que esto se hace sin división o ramificación es bastante inteligente. ¡Pero no hay ningún comentario!Las versiones anteriores de 32 bits de GCC / glibc usaban la
fsin
instrucción, que es sorprendentemente inexacta para algunas entradas. Hay una publicación de blog fascinante que ilustra esto con solo 2 líneas de código .La implementación de fdlibm
sin
en C puro es mucho más simple que la de glibc y está muy bien comentada. Código fuente: fdlibm / s_sin.c y fdlibm / k_sin.cfuente
sin()
; escribagdb a.out
, luegobreak sin
, luegorun
, entoncesdisassemble
.__kernel_sin
está definido en k_sin.c, sin embargo, y es puro C. Haga clic de nuevo; fallé la URL la primera vez.Funciones como seno y coseno se implementan en microcódigo dentro de microprocesadores. Los chips Intel, por ejemplo, tienen instrucciones de montaje para estos. El compilador de CA generará código que llama a estas instrucciones de ensamblaje. (Por el contrario, un compilador de Java no lo hará. Java evalúa las funciones trigonométricas en el software en lugar del hardware, por lo que funciona mucho más lento).
Los chips no usan la serie Taylor para calcular funciones trigonométricas, al menos no del todo. En primer lugar, usan CORDIC , pero también pueden usar una serie corta de Taylor para pulir el resultado de CORDIC o para casos especiales como calcular seno con alta precisión relativa para ángulos muy pequeños. Para obtener más explicaciones, consulte esta respuesta de StackOverflow .
fuente
OK, niños, es hora de los profesionales ... Esta es una de mis mayores quejas con ingenieros de software sin experiencia. Vienen calculando funciones trascendentales desde cero (usando la serie de Taylor) como si nadie hubiera hecho estos cálculos antes en sus vidas. No es verdad. Este es un problema bien definido y ha sido abordado miles de veces por ingenieros de software y hardware muy inteligentes y tiene una solución bien definida. Básicamente, la mayoría de las funciones trascendentales usan polinomios de Chebyshev para calcularlas. En cuanto a qué polinomios se utilizan depende de las circunstancias. Primero, la biblia sobre este asunto es un libro llamado "Computer Approximations" de Hart y Cheney. En ese libro, puede decidir si tiene un sumador de hardware, multiplicador, divisor, etc., y decidir qué operaciones son más rápidas. por ejemplo, si tuvieras un divisor realmente rápido, la forma más rápida de calcular el seno podría ser P1 (x) / P2 (x) donde P1, P2 son polinomios de Chebyshev. Sin el divisor rápido, podría ser solo P (x), donde P tiene muchos más términos que P1 o P2 ... por lo que sería más lento. Entonces, el primer paso es determinar su hardware y lo que puede hacer. Luego, elige la combinación adecuada de polinomios de Chebyshev (generalmente tiene la forma cos (ax) = aP (x) para el coseno, por ejemplo, nuevamente donde P es un polinomio de Chebyshev). Luego decides qué precisión decimal quieres. por ejemplo, si desea una precisión de 7 dígitos, búsquelo en la tabla correspondiente del libro que mencioné, y le dará (para precisión = 7.33) un número N = 4 y un número polinómico 3502. N es el orden de polinomio (entonces es p4.x ^ 4 + p3.x ^ 3 + p2.x ^ 2 + p1.x + p0), porque N = 4. Luego busca el valor real de p4, p3, p2, p1, valores de p0 en la parte posterior del libro por debajo de 3502 (estarán en coma flotante). Luego implementa su algoritmo en software en la forma: (((p4.x + p3) .x + p2) .x + p1) .x + p0 .... y así es como calcularía el coseno a 7 decimales lugares en ese hardware.
Tenga en cuenta que la mayoría de las implementaciones de hardware de operaciones trascendentales en una FPU generalmente implican un microcódigo y operaciones como esta (depende del hardware). Los polinomios de Chebyshev se usan para la mayoría de los trascendentales, pero no para todos. Por ejemplo, la raíz cuadrada es más rápida al usar una iteración doble del método Newton Raphson usando primero una tabla de búsqueda. Una vez más, ese libro "Aproximaciones por computadora" te dirá eso.
Si planea implementar estas funciones, recomendaría a cualquiera que obtenga una copia de ese libro. Realmente es la biblia para este tipo de algoritmos. Tenga en cuenta que hay grupos de medios alternativos para calcular estos valores, como cordics, etc., pero estos tienden a ser mejores para algoritmos específicos donde solo necesita poca precisión. Para garantizar la precisión cada vez, los polinomios de chebyshev son el camino a seguir. Como dije, problema bien definido. Se ha resuelto durante 50 años ... y así es como se hace.
Ahora, dicho esto, existen técnicas mediante las cuales los polinomios de Chebyshev pueden usarse para obtener un resultado de precisión único con un polinomio de bajo grado (como el ejemplo para el coseno anterior). Luego, hay otras técnicas para interpolar entre valores para aumentar la precisión sin tener que recurrir a un polinomio mucho más grande, como el "Método de tablas precisas de Gal". Esta última técnica es a lo que se refiere la publicación que hace referencia a la literatura de ACM. Pero en última instancia, los polinomios de Chebyshev son los que se utilizan para llegar al 90% del camino.
Disfrutar.
fuente
Para
sin
específicamente, utilizando expansión de Taylor le daría:sin (x): = x - x ^ 3/3! + x ^ 5/5! - x ^ 7/7! + ... (1)
seguiría agregando términos hasta que la diferencia entre ellos sea menor que un nivel de tolerancia aceptado o solo por una cantidad finita de pasos (más rápido, pero menos preciso). Un ejemplo sería algo como:
Nota: (1) funciona debido a la aproximación sin (x) = x para ángulos pequeños. Para ángulos más grandes necesita calcular más y más términos para obtener resultados aceptables. Puede usar un argumento while y continuar con cierta precisión:
fuente
Sí, también hay algoritmos de software para calcular
sin
. Básicamente, el cálculo de este tipo de cosas con una computadora digital generalmente se realiza utilizando métodos numéricos como aproximar la serie de Taylor que representa la función.Los métodos numéricos pueden aproximar las funciones a una cantidad arbitraria de precisión y dado que la cantidad de precisión que tiene en un número flotante es finita, se adaptan bastante bien a estas tareas.
fuente
Use la serie Taylor e intente encontrar la relación entre los términos de la serie para no calcular las cosas una y otra vez.
Aquí hay un ejemplo para cosinus:
usando esto podemos obtener el nuevo término de la suma usando el ya usado (evitamos el factorial y x 2p )
fuente
Es una pregunta compleja. La CPU similar a Intel de la familia x86 tiene una implementación de hardware de la
sin()
función, pero es parte de la FPU x87 y ya no se usa en el modo de 64 bits (donde se usan registros SSE2). En ese modo, se utiliza una implementación de software.Hay varias implementaciones de este tipo por ahí. Uno está en fdlibm y se usa en Java. Hasta donde yo sé, la implementación de glibc contiene partes de fdlibm y otras partes aportadas por IBM.
Las implementaciones de software de funciones trascendentales, como el
sin()
uso típico de aproximaciones por polinomios, a menudo obtenidas de la serie Taylor.fuente
sin
ycos
que son más rápidas que las instrucciones de hardware en la FPU. Las bibliotecas más simples e ingenuas tienden a usar las instruccionesfsin
yfcos
.FSIN
con total precisión. Estaría muy agradecido si me dices los nombres de esas bibliotecas rápidas, es interesante echar un vistazo.sin()
es aproximadamente dos veces más rápida de lo que sefsin
calcula (precisamente porque se realiza con menos precisión). Tenga en cuenta que se sabe que el x87 tiene un poco menos de precisión real que sus 79 bits anunciados.Los polinomios de Chebyshev, como se menciona en otra respuesta, son los polinomios donde la mayor diferencia entre la función y el polinomio es lo más pequeña posible. Ese es un excelente comienzo.
En algunos casos, el error máximo no es lo que le interesa, sino el error relativo máximo. Por ejemplo, para la función seno, el error cerca de x = 0 debería ser mucho menor que para valores mayores; quieres un pequeño error relativo . Por lo tanto, calcularía el polinomio de Chebyshev para sen x / x, y multiplicaría ese polinomio por x.
A continuación, tienes que descubrir cómo evaluar el polinomio. Desea evaluarlo de tal manera que los valores intermedios sean pequeños y, por lo tanto, los errores de redondeo sean pequeños. De lo contrario, los errores de redondeo podrían ser mucho más grandes que los errores en el polinomio. Y con funciones como la función seno, si eres descuidado, entonces es posible que el resultado que calculas para sin x sea mayor que el resultado para sin y incluso cuando x <y. Por lo tanto, es necesario elegir cuidadosamente el orden de cálculo y calcular los límites superiores para el error de redondeo.
Por ejemplo, sin x = x - x ^ 3/6 + x ^ 5/120 - x ^ 7/5040 ... Si calcula ingenuamente sin x = x * (1 - x ^ 2/6 + x ^ 4 / 120 - x ^ 6/5040 ...), después de que la función de paréntesis, está disminuyendo, y será suceder que si y es el siguiente número más grande para x, entonces a veces sen y será más pequeño que sen x. En cambio, calcule sen x = x - x ^ 3 * (1/6 - x ^ 2/120 + x ^ 4/5040 ...) donde esto no puede suceder.
Al calcular los polinomios de Chebyshev, por lo general, necesita redondear los coeficientes para duplicar la precisión, por ejemplo. Pero aunque un polinomio de Chebyshev es óptimo, el polinomio de Chebyshev con coeficientes redondeados a doble precisión no es el polinomio óptimo con coeficientes de doble precisión.
Por ejemplo, para sin (x), donde necesita coeficientes para x, x ^ 3, x ^ 5, x ^ 7, etc., haga lo siguiente: Calcule la mejor aproximación de sin x con un polinomio (ax + bx ^ 3 + cx ^ 5 + dx ^ 7) con mayor precisión que el doble, luego redondee a a doble precisión, dando A. La diferencia entre ay A sería bastante grande. Ahora calcule la mejor aproximación de (sen x - Ax) con un polinomio (bx ^ 3 + cx ^ 5 + dx ^ 7). Obtiene diferentes coeficientes, porque se adaptan a la diferencia entre a y A. Redondee b para duplicar la precisión B. Luego, aproxima (sin x - Ax - Bx ^ 3) con un polinomio cx ^ 5 + dx ^ 7 y así sucesivamente. Obtendrá un polinomio que es casi tan bueno como el polinomio original de Chebyshev, pero mucho mejor que Chebyshev redondeado a doble precisión.
A continuación, debe tener en cuenta los errores de redondeo en la elección del polinomio. Encontró un polinomio con un error mínimo en el polinomio ignorando el error de redondeo, pero desea optimizar el polinomio más el error de redondeo. Una vez que tenga el polinomio de Chebyshev, puede calcular los límites para el error de redondeo. Digamos que f (x) es su función, P (x) es el polinomio y E (x) es el error de redondeo. No quieres optimizar | f (x) - P (x) |, desea optimizar | f (x) - P (x) +/- E (x) |. Obtendrá un polinomio ligeramente diferente que intenta mantener los errores polinomiales bajos donde el error de redondeo es grande, y relaja los errores polinómicos un poco donde el error de redondeo es pequeño.
Todo esto le permitirá redondear fácilmente errores de como máximo 0.55 veces el último bit, donde +, -, *, / tienen errores de redondeo como máximo 0.50 veces el último bit.
fuente
En cuanto a la función trigonométrica como
sin()
,cos()
,tan()
no ha habido ninguna mención, después de 5 años, de un aspecto importante de las funciones trigonométricas de alta calidad: reducción en el alcance .Un primer paso en cualquiera de estas funciones es reducir el ángulo, en radianes, a un intervalo de 2 * π. Pero π es irracional, por lo que las reducciones simples como
x = remainder(x, 2*M_PI)
introducir error comoM_PI
, o máquina pi, es una aproximación de π. Entonces, ¿cómo hacerx = remainder(x, 2*π)
?Las primeras bibliotecas utilizaron una precisión extendida o programación diseñada para dar resultados de calidad, pero aún en un rango limitado de
double
. Cuando se solicitó un valor grande comosin(pow(2,30))
, los resultados no tenían sentido o0.0
tal vez con un indicador de error establecido en algo comoTLOSS
pérdida total de precisión oPLOSS
pérdida parcial de precisión.La buena reducción del rango de valores grandes a un intervalo como -π a π es un problema desafiante que rivaliza con los desafíos de la función trigonométrica básica, como
sin()
, en sí misma.Un buen informe es la reducción de argumentos para grandes argumentos: bueno hasta el último momento (1992). Cubre bien el asunto: se analiza la necesidad y cómo eran las cosas en varias plataformas (SPARC, PC, HP, 30 + otros) y proporciona un algoritmo de solución del da resultados de calidad para todos los
double
de-DBL_MAX
aDBL_MAX
.Si los argumentos originales están en grados, pero pueden ser de gran valor, use
fmod()
primero para mejorar la precisión. Un bienfmod()
no introducirá ningún error y proporcionará una excelente reducción del rango.Diversas identidades trigonométricas y
remquo()
ofrecen aún más mejoras. Muestra: sind ()fuente
La implementación real de las funciones de la biblioteca depende del compilador específico y / o proveedor de la biblioteca. Ya sea que se haga en hardware o software, ya sea una expansión Taylor o no, etc., variará.
Me doy cuenta de que eso no es de ninguna ayuda.
fuente
Por lo general, se implementan en software y en la mayoría de los casos no utilizarán las llamadas de hardware correspondientes (es decir, ensamblaje). Sin embargo, como Jason señaló, estos son específicos de la implementación.
Tenga en cuenta que estas rutinas de software no son parte de las fuentes del compilador, sino que se encontrarán en la biblioteca correspondiente, como clib o glibc para el compilador GNU. Ver http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions
Si desea un mayor control, debe evaluar cuidadosamente lo que necesita exactamente. Algunos de los métodos típicos son la interpolación de tablas de búsqueda, la llamada de ensamblaje (que a menudo es lenta) u otros esquemas de aproximación como Newton-Raphson para raíces cuadradas.
fuente
Si desea una implementación en software, no en hardware, el lugar para buscar una respuesta definitiva a esta pregunta es el Capítulo 5 de Recetas numéricas . Mi copia está en una caja, así que no puedo dar detalles, pero la versión corta (si no recuerdo mal) es que tomas
tan(theta/2)
como tu operación primitiva y calculas a los demás desde allí. El cálculo se realiza con una aproximación en serie, pero es algo que converge mucho más rápidamente que una serie de Taylor.Lo siento, no puedo recordar más sin poner mi mano en el libro.
fuente
No hay nada como golpear la fuente y ver cómo alguien realmente lo ha hecho en una biblioteca de uso común; Veamos una implementación de biblioteca C en particular. Elegí uLibC.
Aquí está la función pecado:
http://git.uclibc.org/uClibc/tree/libm/s_sin.c
que parece que maneja algunos casos especiales, y luego lleva a cabo una reducción de argumento para asignar la entrada al rango [-pi / 4, pi / 4], (dividiendo el argumento en dos partes, una gran parte y una cola) antes de llamar
http://git.uclibc.org/uClibc/tree/libm/k_sin.c
que luego opera en esas dos partes. Si no hay cola, se genera una respuesta aproximada utilizando un polinomio de grado 13. Si hay una cola, se obtiene una pequeña suma correctiva basada en el principio de que
sin(x+y) = sin(x) + sin'(x')y
fuente
Cada vez que se evalúa una función de este tipo, en algún nivel es muy probable que:
Si no hay soporte de hardware, entonces el compilador probablemente usa el último método, emitiendo solo código de ensamblador (sin símbolos de depuración), en lugar de usar la biblioteca ac, lo que hace que sea difícil rastrear el código real en su depurador.
fuente
Como mucha gente señaló, depende de la implementación. Pero hasta donde entiendo su pregunta, estaba interesado en una implementación real de software de funciones matemáticas, pero simplemente no logró encontrar una. Si este es el caso, aquí está:
dosincos.c
ubicado en la carpeta glibc root \ sysdeps \ ieee754 \ dbl-64 desempaquetadaTambién puede echar un vistazo a los archivos con la
.tbl
extensión, su contenido no es más que enormes tablas de valores precalculados de diferentes funciones en forma binaria. Es por eso que la implementación es tan rápida: en lugar de calcular todos los coeficientes de cualquier serie que utilicen, solo realizan una búsqueda rápida, que es mucho más rápida. Por cierto, utilizan series Tailor para calcular seno y coseno.Espero que esto ayude.
fuente
Trataré de responder para el caso de
sin()
un programa C, compilado con el compilador C de GCC en un procesador x86 actual (digamos un Intel Core 2 Duo).En el lenguaje C de la biblioteca C estándar incluye funciones matemáticas comunes, no incluidos en el propio idioma (por ejemplo
pow
,sin
ycos
por el poder, seno, coseno y respectivamente). Los encabezados de los cuales se incluyen en matemáticas.h .Ahora en un sistema GNU / Linux, estas funciones de bibliotecas son proporcionadas por glibc (GNU libc o GNU C Library). Pero el compilador de GCC quiere que se vincule al biblioteca matemática (
libm.so
) utilizando el-lm
indicador del compilador para permitir el uso de estas funciones matemáticas.No estoy seguro de por qué no es parte de la biblioteca estándar de C.Sería una versión de software de las funciones de coma flotante, o "soft-float".Aparte: la razón para tener las funciones matemáticas separadas es histórica, y solo tenía la intención de reducir el tamaño de los programas ejecutables en sistemas Unix muy antiguos, posiblemente antes de que las bibliotecas compartidas estuvieran disponibles, hasta donde yo sé.
Ahora el compilador puede optimizar la función estándar de la biblioteca C
sin()
(proporcionada porlibm.so
) para ser reemplazada con una llamada a una instrucción nativa a la función sin () incorporada de su CPU / FPU, que existe como una instrucción FPU (FSIN
para x86 / x87) en procesadores más nuevos como la serie Core 2 (esto es correcto desde el i486DX). Esto dependería de los indicadores de optimización pasados al compilador gcc. Si se le dijera al compilador que escribiera código que se ejecutaría en cualquier procesador i386 o más nuevo, no haría tal optimización. La-mcpu=486
bandera informaría al compilador que era seguro hacer tal optimización.Ahora, si el programa ejecutó la versión de software de la función sin (), lo haría en función de un algoritmo CORDIC (Computadora digital de rotación coordinada) o BKM , o más probablemente en un cálculo de tabla o serie de potencia que se usa comúnmente ahora para calcular tales funciones trascendentales. [Src: http://en.wikipedia.org/wiki/Cordic#Application]
Cualquier versión reciente (desde 2.9x aprox.) De gcc también ofrece una versión integrada de sin,
__builtin_sin()
que se utilizará para reemplazar la llamada estándar a la versión de la biblioteca C, como una optimización.Estoy seguro de que es tan claro como el barro, pero espero que le brinde más información de la que esperaba, y muchos puntos de partida para aprender más usted mismo.
fuente
Si desea ver la implementación real de GNU de esas funciones en C, consulte la última troncal de glibc. Ver la Biblioteca de C de GNU .
fuente
No uses la serie Taylor. Los polinomios de Chebyshev son más rápidos y más precisos, como lo señalan un par de personas arriba. Aquí hay una implementación (originalmente de la ROM ZX Spectrum): https://albertveli.wordpress.com/2015/01/10/zx-sine/
fuente
Calcular el seno / coseno / tangente es realmente muy fácil de hacer a través del código usando la serie Taylor. Escribir uno usted mismo toma como 5 segundos.
Todo el proceso se puede resumir con esta ecuación aquí:
Aquí hay algunas rutinas que escribí para C:
fuente
Versión mejorada del código de la respuesta de Blindy
fuente
La esencia de cómo lo hace radica en este extracto del Análisis numérico aplicado de Gerald Wheatley:
Algunos puntos para mencionar en lo anterior es que algunos algoritmos se interpolan de una tabla, aunque solo para las primeras iteraciones. También tenga en cuenta cómo menciona que las computadoras utilizan polinomios aproximados sin especificar qué tipo de polinomio aproximado. Como han señalado otros en el hilo, los polinomios de Chebyshev son más eficientes que los polinomios de Taylor en este caso.
fuente
si quieres
sin
entoncessi quieres
cos
entoncessi quieres
sqrt
entoncesEntonces, ¿por qué usar un código inexacto cuando las instrucciones de la máquina funcionarán?
fuente