¿Cómo calcula C el pecado () y otras funciones matemáticas?

He estado estudiando detenidamente los desensamblados de .NET y el código fuente de GCC, pero parece que no puedo encontrar en ninguna parte la implementación real de sin() y otras funciones matemáticas … siempre parecen referirse a otra cosa.

¿Alguien puede ayudarme a encontrarlos? Siento que es poco probable que TODO hardware en el que se ejecute C soporte funciones trigonométricas en el hardware, por lo que debe haber un algoritmo de software en alguna parte , ¿no?


Conozco varias formas 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 acerca de cuán reales son los lenguajes de producción, ya que todas mis implementaciones son siempre de varios órdenes de magnitud más lentas, aunque creo que mis algoritmos son bastante ingeniosos (obviamente no lo son).

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 llamas a sin() en un sistema Linux típico x86-64. Aparentemente es más rápido que las instrucciones de ensamblaje fsin . 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 es preciso en todo el rango de valores x , por lo que la biblioteca implementa muchos algoritmos diferentes y su primer trabajo es observar xy decidir qué algoritmo usar. En algunas regiones usa lo que parece ser la serie familiar de Taylor. Varios de los algoritmos primero calculan un resultado rápido, luego, si eso no es lo suficientemente preciso, descartarlo y recurrir a un algoritmo más lento.

Las versiones más antiguas de 32 bits de GCC / glibc usaban la instrucción fsin , que es sorprendentemente inexacta para algunas entradas. Hay una fascinante publicación de blog que ilustra esto con solo 2 líneas de código .

La implementación de fdlibm de sin en C puro es mucho más simple que glibc y está muy bien comentado. Código fuente: fdlibm / s_sin.c y fdlibm / k_sin.c

OK kiddies, time for the pros …. Esta es una de mis mayores quejas con inexpertos ingenieros de software. 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 Chenyshev Polynomials 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” por Hart y Cheney. En ese libro, puede decidir si tiene un sumdor de hardware, multiplicador, divisor, etc., y decidir qué operaciones son más rápidas. por ej., si tuviera 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 apropiada de polinomios de Chebyshev (generalmente tiene la forma cos (ax) = aP (x) para coseno, por ejemplo, donde P es un polinomio de Chebyshev). Luego, usted decide qué precisión decimal quiere. por ejemplo, si desea 7 dígitos de precisión, busque eso en la tabla correspondiente en el libro que mencioné, y le dará (para precisión = 7.33) un número N = 4 y un número polinomial 3502. N es el orden del polinomio (entonces es p4.x ^ 4 + p3.x ^ 3 + p2.x ^ 2 + p1.x + p0), porque N = 4. Luego, busca el valor real de los valores p4, p3, p2, p1, p0 en la parte posterior del libro en 3502 (estarán en coma flotante). Luego implementa su algoritmo en el software en la forma: (((p4.x + p3) .x + p2) .x + p1) .x + p0 …. y así es como se calculó el coseno a 7 decimal lugares en ese hardware.

Tenga en cuenta que la mayoría de las implementaciones de hardware de operaciones trascendentales en una FPU generalmente implican algunos microcódigos 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 para usar una iteración doble del método Newton Raphson utilizando 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 de este tipo de algoritmos. Tenga en cuenta que hay muchos medios alternativos para calcular estos valores, como cordones, etc., pero estos tienden a ser mejores para algoritmos específicos donde solo necesita baja precisión. Para garantizar la precisión en todo momento, los polinomios chebyshev son el camino a seguir. Como dije, problema bien definido. Ha sido resuelto hace 50 años … y así es como se hace.

Ahora, dicho esto, hay técnicas mediante las cuales los polinomios de Chebyshev pueden usarse para obtener un resultado de precisión simple con un polinomio de bajo grado (como el ejemplo del coseno anterior). Luego, hay otras técnicas para interpolar entre los valores para boost la precisión sin tener que ir a un polinomio mucho más grande, como el “Método de Tablas Exactas 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 obtener el 90% del camino hasta allí.

Disfrutar.

Las funciones como seno y coseno se implementan en microcódigo dentro de microprocesadores. Los chips Intel, por ejemplo, tienen instrucciones de ensamblaje para estos. El comstackdor de CA generará un código que llama a estas instrucciones de ensamblaje. (Por el contrario, un comstackdor de Java no lo hará. Java evalúa las funciones trigonométricas en el software en lugar de en el hardware, por lo que se ejecuta mucho más lento).

Las fichas no usan la serie de Taylor para calcular las funciones trigonométricas, al menos no del todo. En primer lugar, utilizan CORDIC , pero también pueden utilizar una serie corta de Taylor para pulir el resultado de CORDIC o para casos especiales, como el cálculo del seno con una precisión relativamente alta para angularjs muy pequeños. Para obtener más explicaciones, consulte esta respuesta de StackOverflow .

Sí, también hay algoritmos de software para calcular el 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 adapta bastante bien a estas tareas.

Es una pregunta compleja. La CPU similar a Intel de la familia x86 tiene una implementación de hardware de la función sin() , pero es parte de la FPU x87 y ya no se usa en el modo de 64 bits (donde en su lugar se usan los registros SSE2). En ese modo, se utiliza una implementación de software.

Hay varias implementaciones similares 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 sin() normalmente usan aproximaciones por polinomios, a menudo obtenidas de series de Taylor.

usa la serie taylor e intenta encontrar la relación entre los términos de la serie para que no puedas calcular las cosas una y otra vez

aquí hay un ejemplo para cosinus:

 double cosinus(double x,double prec) { double t , s ; int p; p = 0; s = 1.0; t = 1.0; while(fabs(t/s) > prec) { p++; t = (-t * x * x) / ((2 * p - 1) * (2 * p)); s += t; } return s;} 

usando esto podemos obtener el nuevo término de la sum usando el ya usado (evitamos el factorial y x ^ 2p)

explicación http://img514.imageshack.us/img514/1959/82098830.jpg

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 más pequeño que para valores más grandes; quieres un pequeño error relativo Entonces, usted calcularía el polinomio de Chebyshev para sin x / x, y multiplicaría ese polinomio por x.

A continuación, debes 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 calcule para sen x sea mayor que el resultado para sen e incluso cuando x

Por ejemplo, sen 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 …), entonces esa función entre paréntesis está disminuyendo, y 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 su lugar, calcule sin x = x – x ^ 3 * (1/6 – x ^ 2/120 + x ^ 4/5040 …) donde esto no puede suceder.

Cuando se calculan los polinomios de Chebyshev, por lo general, se necesitan redondear los coeficientes a doble precisión, por ejemplo. ¡Pero si bien 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 necesitas coeficientes para x, x ^ 3, x ^ 5, x ^ 7 etc. haces lo siguiente: Calcula la mejor aproximación de sen x con un polinomio (ax + bx ^ 3 + cx ^ 5 + dx ^ 7) con una precisión superior a la doble, luego una a doble precisión, dando A. La diferencia entre ay A sería bastante grande. Ahora calcule la mejor aproximación de (sin x – Ax) con un polinomio (bx ^ 3 + cx ^ 5 + dx ^ 7). Obtienes diferentes coeficientes, porque se adaptan a la diferencia entre a y A. Ronda b a doble precisión B. Luego aproxima (sin x – Ax – Bx ^ 3) con un polinomio cx ^ 5 + dx ^ 7 y así sucesivamente. Obtendrás un polinomio que es casi tan bueno como el polinomio original de Chebyshev, pero mucho mejor que el Chebyshev redondeado con doble precisión.

A continuación, debe tener en cuenta los errores de redondeo en la elección del polinomio. Encontraste un polinomio con un error mínimo en el polinomio que ignora el error de redondeo, pero deseas optimizar el error de redondeo más el polinomio. Una vez que tenga el polinomio de Chebyshev, puede calcular los límites para el error de redondeo. Diga 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 en un nivel bajo donde el error de redondeo es grande, y relaja los errores polinomiales un poco donde el error de redondeo es pequeño.

Todo esto hará que redondee fácilmente errores de a lo sumo 0,55 veces el último bit, donde +, -, *, / tienen errores de redondeo de como máximo 0,50 veces el último bit.

Para el sin específicamente, usar la expansión de Taylor te daría:

sin (x): = x – x ^ 3/3! + x ^ 5/5! – x ^ 7/7! + … (1)

seguirás 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 así como:

 float sin(float x) { float res=0, pow=x, fact=1; for(int i=0; i<5; ++i) { res+=pow/fact; pow*=-1*x*x; fact*=(2*(i+1))*(2*(i+1)+1); } return res; } 

Nota: (1) funciona debido a la aproximación sin (x) = x para angularjs pequeños. Para angularjs 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:

 double sin (double x){ int i = 1; double cur = x; double acc = 1; double fact= 1; double pow = x; while (fabs(acc) > .00000001 && i < 100){ fact *= ((2*i)*(2*i+1)); pow *= -1 * x*x; acc = pow / fact; cur += acc; i++; } return cur; } 

La implementación real de las funciones de la biblioteca depende del comstackdor específico y / o del proveedor de la biblioteca. Ya sea que se haga en hardware o software, ya sea una expansión de Taylor o no, etc., variará.

Me doy cuenta de que eso no es de ninguna ayuda.

En cuanto a la función trigonométrica como sin() , cos() , tan() no se menciona, después de 5 años, un aspecto importante de las funciones trigonométricas de alta calidad: reducción del scope .

Un primer paso en cualquiera de estas funciones es reducir el ángulo, en radianes, a un rango de un intervalo de 2 * π. Pero π es irracional, por lo que las reducciones simples como x = remainder(x, 2*M_PI) introducen el error como M_PI , o machine pi, es una aproximación de π. Entonces, ¿cómo hacer x = remainder(x, 2*π) ?

Las primeras bibliotecas usaban una precisión extendida o una progtwigción elaborada para obtener resultados de calidad, pero aún en un rango limitado de double . Cuando se solicitó un valor grande como sin(pow(2,30)) , los resultados fueron insignificantes o 0.0 y tal vez con un indicador de error establecido en algo así como pérdida total de precisión de PLOSS o pérdida parcial de precisión PLOSS .

Una buena reducción de 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í mismo.

Un buen informe es Reducción de argumento para argumentos grandes: Bueno hasta el último bit (1992). Cubre bien el problema: analiza la necesidad y cómo estaban las cosas en varias plataformas (SPARC, PC, HP, más de 30) y proporciona un algoritmo de solución que ofrece resultados de calidad para todo el contenido double desde -DBL_MAX hasta DBL_MAX .


Si los argumentos originales están en grados, pero pueden ser de gran valor, use fmod() primero para una precisión mejorada. Un buen fmod() no introducirá ningún error y proporcionará una excelente reducción de scope.

 // sin(degrees2radians(x)) sin(degrees2radians(fmod(x, 360.0))); // -360.0 <= fmod(x,360) <= +360.0 

Diversas identidades trigonométricas y remquo() ofrecen incluso más mejoras. Muestra: sind ()

Normalmente se implementan en software y no usarán las llamadas de hardware correspondientes (es decir, llamadas de emergencia) en la mayoría de los casos. 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 fonts del comstackdor, sino que se encontrarán en la biblioteca correspondiente, como clib o glibc para el comstackdor 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.

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 recuerdo bien) es que tomes tan(theta/2) como tu operación primitiva y computas las otras desde allí. El cálculo se realiza con una aproximación en serie, pero es algo que converge mucho más rápido que una serie de Taylor.

Lo siento, no puedo recordar más sin poner mi mano en el libro.

Como muchas personas señalaron, depende de la implementación. Pero por lo que entiendo su pregunta, usted estaba interesado en una implementación real de software de funciones matemáticas, pero simplemente no logró encontrar una. Si este es el caso, entonces aquí estás:

  • Descargue el código fuente glibc de http://ftp.gnu.org/gnu/glibc/
  • Mire el archivo dosincos.c ubicado en la dosincos.c desempaquetada glibc root \ sysdeps \ ieee754 \ dbl-64
  • Del mismo modo, puedes encontrar implementaciones del rest de la biblioteca matemática, solo busca el archivo con el nombre apropiado

También puede echar un vistazo a los archivos con la extensión .tbl , sus contenidos no son 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 utilizan, solo hacen una búsqueda rápida, que es mucho más rápida. Por cierto, sí usan la serie Tailor para calcular el seno y el coseno.

Espero que esto ayude.

Trataré de responder por el caso de sin() en un progtwig C, comstackdo con el comstackdor C de GCC en un procesador x86 actual (digamos un Intel Core 2 Duo).

En el lenguaje C, la biblioteca estándar C incluye funciones matemáticas comunes, no incluidas en el idioma en sí (por ejemplo, pow , sin y cos para potencia, seno y coseno, respectivamente). Los encabezados de los cuales están incluidos en math.h.

Ahora en un sistema GNU / Linux, estas funciones de bibliotecas son provistas por glibc (GNU libc o GNU C Library). Pero el comstackdor GCC quiere que libm.so enlace a la biblioteca matemática ( libm.so ) utilizando el -lm comstackdor -lm para habilitar el uso de estas funciones matemáticas. No estoy seguro de por qué no forma parte de la biblioteca C estándar. Estas serían una versión de software de las funciones de punto flotante, o “flotación suave”.

Aparte: la razón por la que las funciones matemáticas se separan es histórica, y solo pretendía reducir el tamaño de los progtwigs ejecutables en sistemas Unix muy antiguos, posiblemente antes de que estuvieran disponibles las bibliotecas compartidas, hasta donde yo sé.

Ahora el comstackdor puede optimizar la función de biblioteca C estándar sin() (proporcionada por libm.so ) para ser reemplazada con una llamada a una instrucción nativa a la función sin() libm.so de la 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 la versión i486DX). Esto dependería de los indicadores de optimización pasados ​​al comstackdor gcc. Si se le dijo al comstackdor que escribiera código que se ejecutaría en cualquier procesador i386 o más nuevo, no haría tal optimización. El -mcpu=486 informaría al comstackdor que era seguro realizar dicha optimización.

Ahora bien, si el progtwig ejecutó la versión de software de la función sin (), lo haría basándose en un CORDIC (Computador DIgital de rotación coordinada) o algoritmo BKM , o más probablemente 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%5D

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 con suerte te da más información de la que esperabas, y muchos puntos de partida para aprender más.

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 la 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 unos pocos 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 usando un polinomio de grado 13. Si hay una cola, se obtiene una pequeña adición correctiva basada en el principio de que sin(x+y) = sin(x) + sin'(x')y

Cuando se evalúa dicha función, entonces en algún nivel es más probable que:

  • Una tabla de valores que se interpola (para aplicaciones rápidas e inexactas, por ejemplo, gráficos de computadora)
  • La evaluación de una serie que converge al valor deseado — probablemente no sea una serie taylor, más probablemente algo basado en una cuadratura elegante como Clenshaw-Curtis.

Si no hay soporte de hardware, entonces el comstackdor probablemente use el último método, emitiendo solo código de ensamblador (sin símbolos de depuración), en lugar de usar la biblioteca de CA, lo que dificulta que rastree el código real en su depurador.

Si desea ver la implementación real de GNU de esas funciones en C, consulte la última línea de glibc. Ver la Biblioteca GNU C

La computación sine / cosine / tangent es en realidad 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í: http://sofes.miximages.com/c/546ecab719ce73dfb34a7496c942972b.png

Aquí hay algunas rutinas que escribí para C:

 double _pow(double a, double b) { double c = 1; for (int i=0; i 

Don’t use Taylor series. Chebyshev polynomials are both faster and more accurate, as pointed out by a couple of people above. Here is an implementation (originally from the ZX Spectrum ROM): https://albertveli.wordpress.com/2015/01/10/zx-sine/

if you want sin then asm volatile (“fsin” : “=t”(vsin) : “0”(xrads)); if you want cos then asm volatile (“fcos” : “=t”(vcos) : “0”(xrads)); if you want sqrt then asm volatile (“fsqrt” : “=t”(vsqrt) : “0”(value)); so why use inaccurate code when the machine instructions will do.