una tipos rango que nocion matematicas matematica lineal funciones funcion elementos dominio c math sine

tipos - que es una funcion matematica



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

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

  • 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 de fantasía como Clenshaw-Curtis.

Si no hay soporte de hardware, es probable que el compilador utilice el último método, que emite solo el código del ensamblador (sin símbolos de depuración), en lugar de usar la biblioteca de CA, por lo que es difícil rastrear el código real en su depurador.

He estado estudiando detenidamente los desensamblajes .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 ... parece que siempre hacen referencia a otra cosa.

¿Alguien puede ayudarme a encontrarlos? Siento que es poco probable que TODO el hardware con el que C se ejecute 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 formas en que se pueden calcular las funciones, y he escrito mis propias rutinas para calcular las funciones usando la serie taylor para la diversión. Tengo curiosidad acerca de cómo lo hacen los lenguajes de producción, ya que todas mis implementaciones son siempre varias órdenes de magnitud más lentas, aunque creo que mis algoritmos son bastante inteligentes (obviamente no lo son).

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

  • Descargue el código fuente de glibc desde http://ftp.gnu.org/gnu/glibc/
  • Mire el archivo dosincos.c ubicado en la dosincos.c glibc root / sysdeps / ieee754 / dbl-64 sin empaquetar
  • 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 extensión .tbl , 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 las series que usan, solo hacen una búsqueda rápida, que es mucho más rápida. Por cierto, usan la serie Tailor para calcular el seno y el coseno.

Espero que esto ayude.

Con respecto a la función trigonométrica como sin() , cos() , tan() no se ha mencionado, después de 5 años, un aspecto importante de las funciones trigonométricas de alta calidad: la reducción de rango .

Un paso temprano 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 máquina pi, es una aproximación de π. Entonces, ¿cómo hacer x = remainder(x, 2*π) ?

Las bibliotecas tempranas utilizaron la precisión extendida o la programación elaborada 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 carecían de sentido o 0.0 y quizás con un indicador de error establecido en algo como TLOSS pérdida total de precisión o PLOSS pérdida parcial de precisión.

Una buena reducción de rango de valores grandes a un intervalo como -π a π es un problema desafiante que compite con los desafíos de la función trigonométrica básica, como el propio sin() .

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

Si los argumentos originales están en grados, pero pueden tener un gran valor, fmod() use fmod() para mejorar la precisión. Un buen fmod() no introducirá ningún error y, por lo tanto, proporcionará una excelente reducción de rango.

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

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

De acuerdo, niños, es hora de que los profesionales ... Esta es una de mis mayores quejas con los ingenieros de software inexpertos. Llegan 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 utilizan 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 "Aproximaciones por computadora" de Hart y Cheney. En ese libro, puede decidir si tiene un sumador, multiplicador, divisor, etc. de hardware, y decidir qué operaciones son las más rápidas. por ejemplo, si tenía 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 ... así 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 (por lo general, tiene la forma cos (ax) = aP (x) para el coseno, por ejemplo, donde P es un polinomio de Chebyshev). Entonces tú decides qué precisión decimal quieres. por ejemplo, si desea una precisión de 7 dígitos, 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 de los polinomio (así que es p4.x ^ 4 + p3.x ^ 3 + p2.x ^ 2 + p1.x + p0), porque N = 4. Luego busca el valor real de los valores de p4, p3, p2, p1, p0 en la parte posterior del libro debajo de 3502 (estarán en punto flotante). Luego implementas tu algoritmo en el software en la forma: (((p4.x + p3) .x + p2) .x + p1) .x + p0 .... y así es como calcularías el coseno con 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 involucran algún 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ápido para usar una iteración doble del método de Newton Raphson utilizando primero una tabla de búsqueda. Una vez más, ese libro "Aproximaciones de computadora" le dirá eso.

Si planea implementar estas funciones, le recomendaría a cualquiera que obtenga una copia de ese libro. Realmente es la biblia de este tipo de algoritmos. Tenga en cuenta que existen muchos medios alternativos para calcular estos valores, como los cordicos, etc., pero estos tienden a ser los mejores para algoritmos específicos en los que solo necesita baja precisión. Para garantizar la precisión en todo momento, los polinomios de Chebyshev son el camino a seguir. Como he dicho, problema bien definido. Se ha solucionado durante 50 años ahora ... y así es como se hace.

Ahora, dicho esto, existen técnicas mediante las cuales los polinomios de Chebyshev se pueden usar 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 valores para aumentar la precisión sin tener que ir 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 se refiere a la literatura de ACM. Pero, en última instancia, los polinomios de Chebyshev son los que se utilizan para obtener el 90% de su camino.

Disfrutar.

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, contribuida por IBM. Desde octubre de 2011, este es el código que realmente se ejecuta cuando se llama a sin() en un sistema Linux x86-64 típico. Aparentemente es más rápido que la instrucción de montaje 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 lo más rápido posible y también es preciso en todo el rango de valores de x , por lo que la biblioteca implementa muchos algoritmos diferentes y su primer trabajo es mirar x y decidir qué algoritmo usar. En algunas regiones usa lo que parece la serie familiar de Taylor. Varios de los algoritmos primero calculan un resultado rápido, luego, si eso no es lo suficientemente preciso, descártelos y recurra a un algoritmo más lento.

Las versiones más antiguas de 32 bits de GCC / glibc utilizaban la instrucción fsin , 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 del sin en C pura es mucho más simple que la de glibc y está bien comentada. Código fuente: fdlibm/s_sin.c y fdlibm/k_sin.c

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 se utilizan 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. Que yo sepa, la implementación de glibc contiene partes de fdlibm y otras partes contribuidas por IBM.

Las implementaciones de software de funciones trascendentales, como sin() suelen usar aproximaciones de polinomios, que a menudo se obtienen de series de Taylor.

Intentaré responder por el caso de sin() en 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, la biblioteca estándar de C incluye funciones matemáticas comunes, no incluidas en el lenguaje en sí (por ejemplo, pow , sin y cos para power, sine y cosine, respectivamente). Los encabezados de los cuales se incluyen en math.h

Ahora en un sistema GNU / Linux, estas funciones de bibliotecas son proporcionadas por glibc (GNU libc o GNU C Library). Pero el compilador GCC quiere que se vincule a la biblioteca matemática ( libm.so ) usando la -lm compilador -lm para habilitar el uso de estas funciones matemáticas. No estoy seguro de por qué no es parte de la biblioteca C estándar. Esta sería una versión de software de las funciones de punto flotante, o "soft-float".

Aparte: el motivo por el cual las funciones matemáticas están separadas es histórico, y su único propósito era reducir el tamaño de los programas ejecutables en sistemas Unix muy antiguos, posiblemente antes de que las bibliotecas compartidas estuvieran disponibles, que yo sepa.

Ahora el compilador puede optimizar la función de biblioteca de C estándar sin() (provista por libm.so ) para que sea reemplazada por una llamada a una instrucción nativa a la función sin() libm.so de su CPU / FPU, que existe como una instrucción de FPU ( FSIN para x86 / x87) en procesadores más nuevos como la serie Core 2 (esto es correcto en gran medida desde el i486DX). Esto dependería de los indicadores de optimización pasados ​​al compilador gcc. Si al compilador se le dijera que escribiera código que se ejecutaría en cualquier i386 o procesador más nuevo, no haría tal optimización. El -mcpu=486 informaría al compilador de que era seguro realizar dicha optimización.

Ahora, si el programa ejecutara la versión de software de la función sin (), lo haría en función de un algoritmo CORDIC (Ordenador Digital de Rotación de Coordenadas) o BKM , o más probablemente una tabla o un cálculo de series de potencias 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 usará 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 te ofrezca más información de la que esperabas, y muchos puntos de partida para que aprendas más.

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

Me doy cuenta de que eso no sirve de nada.

Las funciones como seno y coseno se implementan en microcódigo dentro de microprocesadores. Los chips de Intel, por ejemplo, tienen instrucciones de ensamblaje para estos. El compilador de CA generará un 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 se ejecuta mucho más lento).

Los chips no utilizan la serie de Taylor para calcular funciones trigonométricas, al menos no completamente. En primer lugar, utilizan CORDIC , pero también pueden usar una serie corta de Taylor para pulir el resultado de CORDIC o para casos especiales, como el cálculo del seno con una alta precisión relativa para ángulos muy pequeños. Para más explicación, vea esta respuesta de .

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

A continuación tienes que averiguar 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 es descuidado, puede ser posible que el resultado que calcule para el pecado x sea mayor que el resultado para el pecado y aun 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, 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 ocurrirá que si y es el siguiente número más grande para x, a veces sin y será más pequeño que sin x. En su lugar, calcule sin x = x - x ^ 3 * (1/6 - x ^ 2/120 + x ^ 4/5040 ...) donde esto no pueda suceder.

Al calcular los polinomios de Chebyshev, por lo general es necesario redondear los coeficientes para duplicar la precisión, por ejemplo. ¡Pero mientras 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 una precisión mayor 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). Obtienes diferentes coeficientes, porque se adaptan a la diferencia entre a y A. Redondea b a doble 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 error de redondeo polinomial más. 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 hacia abajo donde el error de redondeo es grande, y relaja los errores polinomiales un poco donde el error de redondeo es pequeño.

Todo esto le permitirá obtener fácilmente errores de redondeo de como máximo 0,55 veces el último bit, donde +, -, *, / tendrá errores de redondeo de como máximo 0,50 veces el último bit.

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 del 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 argumentos para asignar la entrada al rango [-pi / 4, pi / 4], (dividiendo el argumento en dos partes, una parte grande 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

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í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 de tiempo y continuar para una 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; }

Por lo general, se implementan en el software y en la mayoría de los casos no usarán las llamadas correspondientes de hardware (es decir, un ensamblaje). Sin embargo, como Jason señaló, estos son específicos de la implementación.

Tenga en cuenta que estas rutinas de software no forman parte de las fuentes del compilador, sino que se encontrarán en la biblioteca correspondiente como el clib, o glibc para el compilador GNU. Consulte 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 las tablas de consulta, la llamada al ensamblaje (que a menudo es lenta) u otros esquemas de aproximación como Newton-Raphson para raíces cuadradas.

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 mediante métodos numéricos, como la aproximación de la serie de Taylor que representa la función.

Los métodos numéricos pueden aproximar funciones a una cantidad arbitraria de precisión y, como la cantidad de precisión que tiene en un número flotante es finita, se adaptan bastante bien a estas tareas.

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 tomas tan(theta/2) como tu operación primitiva y computas a los demás desde allí. El cálculo se realiza con una serie aproximada, 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.

Si desea ver la implementación GNU real de esas funciones en C, consulte la última troncal de glibc. Vea la biblioteca de GNU C.

use la serie taylor e intente encontrar una relación entre los términos de la serie para que no calcule 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 http://img514.imageshack.us/img514/1959/82098830.jpg

Calcular el seno / coseno / tangente es realmente muy fácil de hacer a través del código utilizando la serie de Taylor. Escribir uno mismo toma como 5 segundos.

El proceso completo se puede resumir con esta ecuación aquí: http://upload.wikimedia.org/math/5/4/6/546ecab719ce73dfb34a7496c942972b.png

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)); }

si quieres el pecado entonces asm volatile ("fsin": "= t" (vsin): "0" (xrads)); si desea cos, asm volatile ("fcos": "= t" (vcos): "0" (xrads)); si desea sqrt, asm volatile ("fsqrt": "= t" (vsqrt): "0" (valor)); Entonces, ¿por qué usar un código inexacto cuando las instrucciones de la máquina funcionarán?