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

248

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).

Madeja
fuente
2
Tenga en cuenta que esta implementación depende. Debe especificar en qué implementación está más interesado.
jason
3
Etiqueté .NET y C porque busqué en ambos lugares y tampoco pude entenderlo. Aunque mirando el desmontaje de .NET parece que podría estar llamando a C no administrado, por lo que yo sé, tienen la misma implementación.
Hank

Respuestas:

213

En GNU libm, la implementación de sindepende 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 las fsininstrucciones 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) == xes 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 sino cosno.

  • 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.

double t = (x * hpinv + toint);
double xn = t - toint;

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 fsininstrucció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 sinen 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.c

Jason Orendorff
fuente
35
Para ver que este es realmente el código que se ejecuta en x86: compile un programa que llame sin(); escriba gdb a.out, luego break sin, luego run, entonces disassemble.
Jason Orendorff
55
@ Henry: no cometas el error de pensar que es un buen código. Es realmente terrible , ¡no aprendas a codificar de esa manera!
Thomas Bonini
2
@Andreas Hmm, tienes razón, el código de IBM se ve bastante horrible en comparación con fdlibm. Edité la respuesta para agregar enlaces a la rutina sinusoidal de fdlibm.
Jason Orendorff
3
@Henry: __kernel_sinestá definido en k_sin.c, sin embargo, y es puro C. Haga clic de nuevo; fallé la URL la primera vez.
Jason Orendorff
3
El código sysdeps vinculado es particularmente interesante porque está correctamente redondeado. Es decir, aparentemente da la mejor respuesta posible para todos los valores de entrada, que solo ha sido posible recientemente. En algunos casos, esto puede ser lento porque puede ser necesario calcular muchos dígitos adicionales para garantizar el redondeo correcto. En otros casos, es extremadamente rápido: para números lo suficientemente pequeños, la respuesta es solo el ángulo.
Bruce Dawson
66

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 .

John D. Cook
fuente
10
Las funciones matemáticas trascendentales como seno y coseno pueden implementarse en microcódigo o como instrucciones de hardware en procesadores de escritorio y servidor actuales de 32 bits. Este no fue siempre el caso, hasta el i486 (DX) todos los cálculos de coma flotante se realizaron en software ("soft-float") para la serie x86 sin un coprocesador separado. No todos los cuales (FPU) incluían funciones trascendentales (por ejemplo, Weitek 3167).
mctylr
1
¿Puedes ser mas específico? ¿Cómo se "pule" una aproximación usando una serie de Taylor?
Hank
44
En cuanto a "pulir" una respuesta, suponga que está calculando seno y coseno. Suponga que conoce el valor exacto de ambos en un punto (por ejemplo, de CORDIC) pero desea el valor en un punto cercano. Luego, para una pequeña diferencia h, puede aplicar las aproximaciones de Taylor f (x + h) = f (x) + h f '(x) o f (x + h) = f (x) + h f' (x) + h ^ 2 f '' (x) / 2.
John D. Cook
66
Los chips x86 / x64 tienen una instrucción de ensamblaje para calcular el seno (fsin), pero esta instrucción a veces es bastante inexacta y, por lo tanto, rara vez se usa más. Ver randomascii.wordpress.com/2014/10/09/… para más detalles. La mayoría de los otros procesadores no tienen instrucciones para el seno y el coseno porque calcularlos en el software proporciona más flexibilidad e incluso puede ser más rápido.
Bruce Dawson
3
Las cosas cordicas dentro de los chips de inteligencia generalmente NO se usan. Primero, la precisión y resolución de la operación es extremadamente importante para muchas aplicaciones. Cordic es notoriamente inexacto cuando llegas al séptimo dígito más o menos, e impredecible. En segundo lugar, escuché que hay un error en su implementación, lo que causa aún más problemas. Eché un vistazo a la función sin para linux gcc y, efectivamente, usa chebyshev. No se utiliza el material incorporado. Ah, también, el algoritmo cordic en el chip es más lento que la solución de software.
Donald Murray
63

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.

Donald Murray
fuente
66
No podría estar más de acuerdo con las primeras oraciones. Además, vale la pena recordar que calcular funciones especiales con precisión garantizada es un problema difícil . Las personas inteligentes que mencionas pasan la mayor parte de su vida haciendo esto. Además, en una nota más técnica, los polinomios min-max son los buscados, y los polinomios de Chebyshev son representantes más simples para ellos.
Alexandre C.
161
-1 para el tono poco profesional y laberíntica (y suavemente groseros), y por el hecho de que el no redundante real contenido de esta respuesta, despojado de la laberíntica y condescendencia, básicamente se reduce a "A menudo utilizan polinomios de Chebyshev; ver este libro para más detalles, ¡es realmente bueno! " Lo cual, ya sabes, bien podría ser absolutamente correcto, pero no es realmente el tipo de respuesta autónoma que queremos aquí en SO. Sin embargo, condensado así, habría hecho un comentario decente sobre la pregunta.
Ilmari Karonen
2
En los primeros años de desarrollo del juego, generalmente se hacía con tablas de búsqueda (necesidad crítica de velocidad). Por lo general, no usamos las funciones lib estándar para esas cosas.
topspin
44
Utilizo tablas de búsqueda en sistemas integrados con bastante frecuencia y bitianos (en lugar de radianes), pero esto es para una aplicación especializada (como sus juegos). Creo que el tipo está interesado en cómo el compilador de c calcula el pecado para números de coma flotante ...
Donald Murray
1
Ah, hace 50 años. Empecé a jugar con eso en la serie Burroughs B220 con McLaren. Más tarde hardware CDC y luego Motorola 68000. Arcsin estaba desordenado: elegí el cociente de dos polinomios y desarrollé código para encontrar los coeficientes óptimos.
Rick James
15

Para sinespecí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:

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 á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:

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;

}
Ciego
fuente
1
Si modifica un poco los coeficientes (y los codifica en un polinomio), puede detener aproximadamente 2 iteraciones antes.
Rick James
14

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.

Mehrdad Afshari
fuente
12
Una implementación real probablemente no usará una serie Taylor, ya que hay formas más eficientes. Solo necesita aproximarse correctamente en el dominio [0 ... pi / 2], y hay funciones que ofrecerán una buena aproximación de manera más eficiente que una serie de Taylor.
David Thornley
2
@David: estoy de acuerdo. Fui lo suficientemente cuidadoso como para mencionar la palabra "me gusta" en mi respuesta. Pero la expansión de Taylor es simple para explicar la idea detrás de los métodos que se aproximan a las funciones. Dicho esto, he visto implementaciones de software (no estoy seguro si fueron optimizadas) que usaban la serie Taylor.
Mehrdad Afshari
1
En realidad, las aproximaciones polinómicas son una de las formas más eficientes de calcular funciones trigonométricas.
Jeremy Salwen
13

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:

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 suma usando el ya usado (evitamos el factorial y x 2p )

explicación

Hannoun Yassir
fuente
2
¿Sabía que puede usar la API de Google Chart para crear fórmulas como esta usando TeX? code.google.com/apis/chart/docs/gallery/formulas.html
Gab Royer
11

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.

Thomas Pornin
fuente
3
Los registros SSE2 no se utilizan para calcular sin (), ni en modo x86 ni en modo x64 y, por supuesto, sin se calcula en hardware independientemente del modo. Hola, es 2010 vivimos en :)
Igor Korkhov
77
@ Igor: eso depende de qué biblioteca matemática estés viendo. Resulta que las bibliotecas de matemáticas más optimizadas en x86 usan implementaciones de software SSE siny cosque son más rápidas que las instrucciones de hardware en la FPU. Las bibliotecas más simples e ingenuas tienden a usar las instrucciones fsiny fcos.
Stephen Canon
@Stephen Canon: ¿tienen esas bibliotecas rápidas una precisión de 80 bits como lo hacen los registros de FPU? Tengo una sospecha muy disimulada de que favorecen la velocidad sobre la precisión, lo que por supuesto es razonable en muchos escenarios, por ejemplo en los juegos. Y sí creo que calcular el seno con una precisión de 32 bits usando SSE y tablas intermedias precalculadas podría ser más rápido que si se usa FSINcon total precisión. Estaría muy agradecido si me dices los nombres de esas bibliotecas rápidas, es interesante echar un vistazo.
Igor Korkhov
@Igor: en x86 en modo de 64 bits, al menos en todos los sistemas tipo Unix que conozco, la precisión está limitada a 64 bits, no a los 79 bits de la FPU x87. La implementación del software sin()es aproximadamente dos veces más rápida de lo que se fsincalcula (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.
Thomas Pornin
1
De hecho, las implementaciones de sin () de 32 bits y 64 bits en las bibliotecas de tiempo de ejecución msvc no utilizan la instrucción FSIN. De hecho, dan resultados diferentes, por ejemplo, sin (0.70444454416678126). Esto dará como resultado 0.64761068800896837 (con una tolerancia de 0.5 * (eps / 2)) en un programa de 32 bits, y dará como resultado 0.64761068800896848 (incorrecto) en uno de 64 bits.
e.tadeu
9

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.

gnasher729
fuente
1
Esta es una buena explicación de cómo se puede calcular sin (x) de manera eficiente, pero realmente no parece responder a la pregunta del OP, que es específicamente sobre cómo las bibliotecas / compiladores de C comunes lo calculan.
Ilmari Karonen
Los polinomios de Chebyshev minimizan el valor absoluto máximo durante un intervalo, pero no minimizan la mayor diferencia entre una función objetivo y el polinomio. Los polinomios de Minimax hacen eso.
Eric Postpischil
9

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 como M_PI, o máquina pi, es una aproximación de π. Entonces, ¿cómo hacer x = 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 como sin(pow(2,30)), los resultados no tenían sentido o 0.0tal vez con un indicador de error establecido en algo como TLOSSpérdida total de precisión o PLOSSpé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_MAXa DBL_MAX.


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

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

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

chux - Restablece a Monica
fuente
6

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.

John Bode
fuente
5

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.

mnemosyn
fuente
5

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.

Norman Ramsey
fuente
5

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 quesin(x+y) = sin(x) + sin'(x')y

Moschops
fuente
4

Cada vez que se evalúa una función de este tipo, en algún nivel es muy probable que:

  • Una tabla de valores que se interpola (para aplicaciones rápidas e imprecisas, por ejemplo, gráficos por 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 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.

James
fuente
4

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á:

  • Descargue el código fuente de glibc de http://ftp.gnu.org/gnu/glibc/
  • Mire el archivo dosincos.cubicado en la carpeta glibc root \ sysdeps \ ieee754 \ dbl-64 desempaquetada
  • De manera similar, puede encontrar implementaciones del resto de la biblioteca matemática, solo busque el archivo con el nombre apropiado

También puede echar un vistazo a los archivos con la .tblextensió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.

Igor Korkhov
fuente
4

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, siny cospor 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 -lmindicador 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 por libm.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 ( FSINpara 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=486bandera 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.

mctylr
fuente
3

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 .

Chris Tonkinson
fuente
3

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/

Albert Veli
fuente
2
Esto realmente no parece responder a la pregunta que se le hizo. El OP pregunta cómo se calculan las funciones trigonométricas mediante compiladores / bibliotecas comunes de C (y estoy bastante seguro de que ZX Spectrum no califica), no cómo deberían calcularse. Sin embargo, este podría haber sido un comentario útil sobre algunas de las respuestas anteriores.
Ilmari Karonen
1
Ah, tienes razon. Debería haber sido un comentario y no una respuesta. No he usado SO en mucho tiempo y olvidé cómo funciona el sistema. De todos modos, creo que la implementación de Spectrum es relevante porque tenía una CPU realmente lenta y la velocidad era esencial. Entonces, el mejor algoritmo sigue siendo bastante bueno, por lo que sería una buena idea que las bibliotecas C implementen funciones trigonométricas utilizando polinomios de Chebyshev.
Albert Veli
2

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í:

expansión del pecado y el costo

Aquí hay algunas rutinas que escribí para C:

double _pow(double a, double b) {
    double c = 1;
    for (int i=0; i<b; i++)
        c *= a;
    return c;
}

double _fact(double x) {
    double ret = 1;
    for (int i=1; i<=x; i++) 
        ret *= i;
    return ret;
}

double _sin(double x) {
    double y = x;
    double s = -1;
    for (int i=3; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _cos(double x) {
    double y = 1;
    double s = -1;
    for (int i=2; i<=100; i+=2) {
        y+=s*(_pow(x,i)/_fact(i));
        s *= -1;
    }  
    return y;
}
double _tan(double x) {
     return (_sin(x)/_cos(x));  
}
usuario1432532
fuente
44
Esta es una implementación bastante mala ya que no utiliza que los términos sucesivos de las series seno y coseno tengan cocientes muy simples. Lo que significa que uno puede reducir el número de multiplicaciones y divisiones de O (n ^ 2) aquí a O (n). Se logran reducciones adicionales dividiendo a la mitad y cuadrando, por ejemplo, como se hace en la biblioteca matemática bc (calculadora de multiprecisión POSIX).
Lutz Lehmann
2
Tampoco parece estar respondiendo la pregunta como se le preguntó; El OP pregunta cómo se calculan las funciones trigonométricas por los compiladores / bibliotecas comunes de C, no para las implementaciones personalizadas.
Ilmari Karonen
2
Creo que es una buena respuesta, ya que responde al espíritu de la pregunta que (y solo puedo adivinar, por supuesto) la curiosidad sobre una función de "caja negra" como el pecado (). Es la única respuesta aquí que le da a uno la oportunidad de comprender rápidamente lo que está sucediendo al pasar por alto en unos segundos en lugar de leer un código fuente C optimizado.
Mike M
de hecho, las bibliotecas usan la versión mucho más optimizada, al darse cuenta de que una vez que tiene un término, puede obtener el siguiente término multiplicando algunos valores. Vea un ejemplo en la respuesta de Blindy . Estás calculando la potencia y los factoriales una y otra vez, lo que es mucho más lento
phuclv
0

Versión mejorada del código de la respuesta de Blindy

#define EPSILON .0000000000001
// this is smallest effective threshold, at least on my OS (WSL ubuntu 18)
// possibly because factorial part turns 0 at some point
// and it happens faster then series element turns 0;
// validation was made against sin() from <math.h>
double ft_sin(double x)
{
    int k = 2;
    double r = x;
    double acc = 1;
    double den = 1;
    double num = x;

//  precision drops rapidly when x is not close to 0
//  so move x to 0 as close as possible
    while (x > PI)
        x -= PI;
    while (x < -PI)
        x += PI;
    if (x > PI / 2)
        return (ft_sin(PI - x));
    if (x < -PI / 2)
        return (ft_sin(-PI - x));
//  not using fabs for performance reasons
    while (acc > EPSILON || acc < -EPSILON)
    {
        num *= -x * x;
        den *= k * (k + 1);
        acc = num / den;
        r += acc;
        k += 2;
    }
    return (r);
}
Rugnar
fuente
0

La esencia de cómo lo hace radica en este extracto del Análisis numérico aplicado de Gerald Wheatley:

Cuando su programa de software le pide a la computadora que obtenga un valor de ingrese la descripción de la imagen aquío ingrese la descripción de la imagen aquí, ¿se ha preguntado cómo puede obtener los valores si las funciones más poderosas que puede calcular son polinomios? ¡No busca estos en tablas e interpola! Por el contrario, la computadora aproxima cada función que no sea polinomios de algún polinomio que está diseñada para dar los valores con mucha precisión.

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.

Dean P
fuente
-1

si quieres sinentonces

 __asm__ __volatile__("fsin" : "=t"(vsin) : "0"(xrads));

si quieres cosentonces

 __asm__ __volatile__("fcos" : "=t"(vcos) : "0"(xrads));

si quieres sqrtentonces

 __asm__ __volatile__("fsqrt" : "=t"(vsqrt) : "0"(value));

Entonces, ¿por qué usar un código inexacto cuando las instrucciones de la máquina funcionarán?

usuario80998
fuente