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

Bien, chicos, 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 utilizan Chebyshev Polinomios para calcularlos. En cuanto a qué polinomios se utilizan depende de las circunstancias. Primero, la biblia sobre este asunto es un libro llamado "Aproximaciones Computacionales" por Hart y Cheney. En ese libro, puede decidir si tiene un adder de hardware, multiplicador, divisor, etc., y decidir qué operaciones son más rápidas. por ejemplo, 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 mucho más términos que P1 o P2....so sería más lento. Por lo tanto, el primer paso es determinar su hardware y lo que puede hacer. Entonces usted elige la combinación apropiada de polinomios de Chebyshev (es generalmente de la forma cos ( ax) = aP(x) para el coseno por ejemplo, otra vez 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, busque eso en la tabla apropiada en el 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 del polinomio (por lo 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 p4,p3,p2,p1,p0 en la parte posterior del libro bajo 3502 (estarán en coma flotante). Luego implementa su algoritmo en el software en el formulario: (((p4.x + p3).x + p2).x + p1).x + p0 ....y así es como calcularías el coseno a 7 decimales en ese hardware.

Tenga en cuenta que la mayoría del hardware las implementaciones 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ápida para usar una doble iteración del método de Newton raphson usando primero una tabla de búsqueda. Una vez más, ese libro "Aproximaciones de computadora" te dirá eso.

Si planea implementar estas funciones, le recomendaría a cualquiera que obtenga una copia de ese libro. Se realmente es la biblia para este tipo de algoritmos. Tenga en cuenta que hay racimos de medios alternativos para calcular estos valores como cordics, etc, pero estos tienden a ser mejores para algoritmos específicos donde solo se necesita baja 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 desde hace 50 años.....y así es como se hace.

Ahora bien, dicho esto, hay técnicas mediante las cuales el Chebyshev los polinomios se pueden usar para obtener un solo resultado de precisión 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 ir a un polinomio mucho más grande, como el "Método de Tablas Precisas de Gal". Esta última técnica es a la que se refiere el post referido a la literatura de la MCA. Pero en última instancia, los polinomios de Chebyshev son los que se utilizan para obtener el 90% del camino alli.

Disfruta.

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

Los chips no usan series de Taylor para calcular funciones trigonométricas, al menos no del todo. 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 de seno con alta precisión relativa para ángulos muy pequeños. Para más explicación, vea esta respuesta de StackOverflow.

Sí, también hay algoritmos de software para calcular sin. Básicamente, calcular este tipo de cosas con una computadora digital generalmente se hace usando métodos numéricos como aproximar 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 dado que la cantidad de precisión que tiene en un número flotante es finita, se adaptan bastante bien a estas tareas.

Es una cuestión 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 modo de 64 bits (donde se usan registros SSE2 en su lugar). En ese modo, se utiliza una implementación de software.

Existen varias implementaciones de allí. Uno está en fdlibm y se usa en Java. Hasta donde 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() típicamente usan aproximaciones por polinomios, a menudo obtenidas de series de Taylor.

Usa la serie de taylor y trata de 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 cosino:

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 utilizado (evitamos el factorial y x^2p)

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

Los polinomios de Chebyshev, como se mencionó en otra respuesta, son los polinomios donde la mayor diferencia entre la función y el polinomio es lo más pequeña posible. Es un excelente comienzo.

En algunos casos, el error máximo no es lo que le interesa, sino el error máximo relativo. Por ejemplo, para la función seno, el error cerca de x = 0 debe ser mucho menor que para valores más grandes; desea un pequeño error relativo. Así que usted calcularía el Chebyshev polinomio para el pecado x / x, y multiplicar 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 llegar a ser mucho más grandes que los errores en el polinomio. Y con funciones como la función seno, si eres descuidado, entonces puede ser posible que el resultado que calculas para sin x sea mayor que el resultado para el pecado y incluso cuando x

Por ejemplo, sin x = x - x^3/6 + x^5 / 120 - x^7 / 5040... Si se calcula ingenuamente sen 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 mayor a x, entonces a veces el pecado y será menor que el pecado x. En su lugar, calcule el pecado x = x - x^3 * (1/6-x^2 / 120 + x^4/5040...) donde esto no puede suceder.

Al calcular polinomios de Chebyshev, generalmente necesita redondear los coeficientes para duplicar la precisión, por ejemplo. Pero mientras que 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. usted hace lo siguiente: Calcular el mejor aproximación de sin x con un polinomio (ax + bx^3 + cx^5 + dx^7) con mayor precisión que el doble, luego redondear a a doble precisión, dando A. La diferencia entre a y A sería bastante grande. Ahora calcular la mejor aproximación de (sin x-Ax) con un polinomio (b x^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. Entonces aproximado (sin x - Ax - Bx^3) con un polinomio cx^5 + dx^7 y así sucesivamente. Usted obtener 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. Ha encontrado 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. Decir f (x) es su función, P (x) es el polinomio, y E (x) es el error de redondeo. Usted no quiere optimizar | f (x) - P (x) |, desea optimizar | f (x) - P (x) +/- E (x) |. Obtendrá un polinomio ligeramente diferente que intenta mantener los errores polinómicos abajo 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 de como máximo 0.50 veces la última parte.

Para 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 para 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 que necesita para 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;

}

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

En cuanto 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: Reducción de rango.

Un primer paso en cualquiera de estas funciones es reducir el ángulo, en radianes, a un intervalo de 2*π. Pero π es irracional reducciones tan 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 utilizadas extendidas precisión o programación hecha a mano para dar resultados de calidad, pero todavía sobre un rango limitado de double. Cuando se solicitó un valor grande como sin(pow(2,30)), los resultados no tenían sentido o 0.0 y tal vez con una bandera de error establecida en algo como TLOSS pérdida total de precisión o PLOSS pérdida parcial de precisión.

La 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(), sí mismo.

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

Si los argumentos originales están en grados, pero pueden ser de un valor grande, use fmod() primero para mejorar la precisión. Una buena fmod() voluntad introduce sin error y así proporciona una excelente reducción del 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()

Normalmente se implementan en software y no usarán las llamadas de hardware correspondientes (es decir, aseembly) en la mayoría de los casos. Sin embargo, como Jason señaló, estos son específicos de 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 el clib, o glibc para el compilador GNU. Véase http://www.gnu.org/software/libc/manual/html_mono/libc.html#Trig-Functions

Si quieres mayor control, usted debe evaluar cuidadosamente lo que necesita exactamente. Algunos de los métodos típicos son la interpolación de tablas de consulta, la llamada de ensamblado (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 esto) es que tomas tan(theta/2) como tu operación primitiva y calculas las otras a partir de ahí. El cálculo se hace con una aproximación de serie, pero es algo que converge mucho más rápidamente que un Taylor serie.

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 un verdadero software implementación de funciones matemáticas, pero simplemente no logró encontrar una. Si este es el caso, entonces aquí estás:

• Descargue el código fuente de glibc desde http://ftp.gnu.org/gnu/glibc /
• Mira el archivo dosincos.c ubicado en raíz glibc descomprimida \sysdeps\ieee754\carpeta dbl-64
• Del mismo modo se puede encontrar implementaciones del resto de la biblioteca de matemáticas, simplemente 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 cualquier serie que utilicen, simplemente hacen una búsqueda rápida, que es mucho más rápida. Por cierto, usan Sastre serie para calcular seno y coseno.

Espero que esto ayude.

Trataré de 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 C estándar incluye funciones matemáticas comunes, no incluidas en el lenguaje en sí (p. ej. pow, sin y cos para la potencia, seno y coseno respectivamente). Cuyas cabeceras están incluidas en math.h .

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

Aparte: La razón para tener las funciones matemáticas separadas es histórica, y fue simplemente la intención de reducir el tamaño de los programas ejecutables en muy sistemas Unix 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 (FSIN para x86/x87) en procesadores más nuevos como la serie Core 2 (esto es correcto casi desde el i486DX). Esto dependería de los indicadores de optimización pasado al compilador gcc. Si se le dijera al compilador que escribiera código que se ejecutaría en cualquier procesador i386 o más reciente, no haría tal optimización. El indicador -mcpu=486 informaría al compilador que era seguro hacer tal optimización.

Ahora, si el programa ejecutara la versión de software de la función sin (), lo haría basado en un CORDIC (Computadora digital de Rotación de coordenadas) o algoritmo BKM, o más probablemente una tabla o serie de potencias cálculo que se utiliza comúnmente ahora para calcular tales funciones trascendentales. [Src: http://en.wikipedia.org/wiki/Cordic#Application]

Cualquier reciente (desde 2.9 x aprox.) version of gcc also offers a built-in version of sin, __builtin_sin() that it will used to replace the standard call to the C library version, as an optimization.

Estoy seguro de que es tan claro como el barro, pero con suerte te da más información de la que esperabas, y un montón de puntos de salto a aprende más tú mismo.

No hay nada como golpear el código fuente y ver cómo alguien realmente lo ha hecho en una biblioteca de uso común; echemos un vistazo a una implementación de biblioteca de C en particular. Elegí uLibC.

Aquí está la función sin:

Http://git.uclibc.org/uClibc/tree/libm/s_sin.c

Que parece que maneja algunos casos especiales, y luego lleva a cabo alguna reducción de argumentos para mapear 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

Cada vez que se evalúa una función de este tipo, 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 por computadora)
• La evaluación de una serie que converge al valor deseado --- probablemente no una serie de taylor, más probablemente algo basado en una cuadratura de fantasía como Clenshaw-Curtis.

Si no hay soporte de hardware, el compilador probablemente usa este último método, emitiendo solo código de ensamblador (sin símbolos de depuración), en lugar de usar una biblioteca de c, lo que le dificulta rastrear el código real en su depurador.

Si quieres ver la implementación actual de GNU de esas funciones en C, echa un vistazo al último tronco de glibc. Vea la Biblioteca de GNU C.

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

Todo El proceso 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));  
}

No use la serie de Taylor. Los polinomios de Chebyshev son más rápidos y más precisos, como señalaron un par de personas anteriormente. Aquí hay una implementación (originalmente de la ROM de ZX Spectrum): https://albertveli.wordpress.com/2015/01/10/zx-sine /

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