Aproximación de funciones matemáticas

[Picuino]

Los compiladores de c suelen venir con un conjunto de funciones matemáticas estandar para realizar cáculos con ellas.

Las funciones más usuales son las trigonométricas (seno, coseno, tangente, arcoseno, arcocoseno, arcotangente) y las exponenciales y logarítmicas (logaritmo natural, logaritmo base 10, exponencial, 10^x, etc)

Las funciones de la librería estandar vienen con una precisión predefinida, que suele ser la mayor precisión que se puede conseguir con el sistema de numeración utilizado. Por ejemplo, para números en coma flotante de 32bits la precisión máxima es aproximadamente de 7 decimales. En el caso de utilizar números en coma flotante de 64 bits, se puede conseguir una precisión de 16 decimales.

Todas estas funciones presentan ciertos problemas:

1.- Rapidez:
En ocasiones se necesita realizar un cálculo de una función matemática (seno, coseno, etc) en un microcontrolador con mucha rapidez y relativamente baja precisión. Las funciones de la librería estandar son precisas a costa de ser más lentas.
En este caso se desea conseguir una aproximación a la función diferente a la librería estandar: más rápida a costa de la precisión.

2.- Variedad:
En otras ocasiones necesitamos que el microcontrolador calcule funciones matemáticas conocidas, que no están en la librería estandar (por ejemplo funciones estadísticas). En este caso necesitamos que el microcontrolador calcule por sí mismo la función matemática.

3.- Funciones desconocidas:
En ocasiones necesitamos almacenar una función matemática de la cual no tenemos la fórmula. Una solución puede ser utilizar tablas de datos e interpolar. Este sistema tiene el inconveniente de ocupar mucha memoria. Puede llegar a ocupar más memoria de la disponible en el microcontrolador. En este caso se necesita convertir la tabla de datos en una función matemática que la aproxime y que se calcule ocupando menos memoria.


En este hilo voy a presentar algunos recursos para solucionar estos problemas.

Un saludo.

[Picuino]

BIBLIOGRAFÍA: APROXIMACIÓN DE FUNCIONES MATEMÁTICAS POR COMPUTADOR


Abramowitz and Stegun: Handbook of Mathematical Functions
http://people.math.sfu.ca/~cbm/aands/toc.htm
http://www.nr.com/aands/

John F. Hart and Cheney: Computer Approximations. ed John Wiley & Sons
http://books.google.es/books/about/Computer_Approximations.html?id=jPVQAAAAMAAJ&redir_esc=y


StatLib---Applied Statistics algorithms
http://lib.stat.cmu.edu/apstat/
Source codes: http://people.sc.fsu.edu/~jburkardt/c_src/c_src.html

Saludos.

[Picuino]

Ejemplo: Función sen(x) aproximada con cuatro decimales
Tomado de página 76 de A&S

Código: [Seleccionar]

   sin(x) = x - 0.16605*x^3 + 0.00761 *x^5
   sin(x) = x * (1 + x^2  * (-0.16605 + x^2 * 0.00761)))

Lo convertimos a una forma más sencilla para programarlo:

Código: [Seleccionar]

   xx = x * x
   sin = 0.00761
   sin = sin * xx - 0.16605
   sin = sin * xx + 1
   sin = sin * x

Esta función es rápida, con precisión suficiente para muchos casos.
Adjunto gráfico con el error cometido.

Error:

Código: [Seleccionar]

x sin(x) p(x) error relativo
0,00000 0,00000 0,00000 0,00000
0,06283 0,06279 0,06279 0,00000
0,12566 0,12533 0,12533 0,00001
0,18850 0,18738 0,18739 0,00002
0,25133 0,24869 0,24870 0,00004
0,31416 0,30902 0,30903 0,00005
0,37699 0,36812 0,36815 0,00008
0,43982 0,42578 0,42582 0,00010
0,50265 0,48175 0,48181 0,00012
0,56549 0,53583 0,53590 0,00014
0,62832 0,58779 0,58788 0,00015
0,69115 0,63742 0,63753 0,00016
0,75398 0,68455 0,68466 0,00017
0,81681 0,72897 0,72909 0,00017
0,87965 0,77051 0,77063 0,00015
0,94248 0,80902 0,80912 0,00013
1,00531 0,84433 0,84441 0,00010
1,06814 0,87631 0,87636 0,00006
1,13097 0,90483 0,90484 0,00002
1,19381 0,92978 0,92974 -0,00003
1,25664 0,95106 0,95097 -0,00009
1,31947 0,96858 0,96846 -0,00013
1,38230 0,98229 0,98213 -0,00016
1,44513 0,99211 0,99195 -0,00016
1,50796 0,99803 0,99791 -0,00012
1,57080 1,00000 1,00000 0,00000


Saludos.

[migsantiago]

Gracias por compartirlo Picuino. 

[AngelGris]

Muy bueno Picuino, como nos tienes acostumbrados 

[Picuino]

Biblioteca en C de las funciones matemáticas más comunes en precisión doble (64 bits):

Cabecera de librería matemática:
cmath.h

Código: C

1. #define HUGE       (1.701411733192644270e38)
2. #define TWO_PI     (0.63661977236758134307553505349)      /*  2/pi     */
3. #define PI_TWO     (1.57079632679489661923132169163)      /*  pi/2     */
4. #define PI         (3.14159265358979323846264338327)      /*  pi       */
5. #define INV_PI     (0.31830988618379067153776752674)      /*  1/pi     */
6. #define LN10       (2.30258509299404568401799145468)      /* ln(10)    */
7. #define INV_LN10   (0.43429448190325182765112891892)      /* 1/log(10) */
8. #define LN2        (0.69314718055994530941723212146)      /* log(2)    */
9. #define INV_LN2    (1.44269504088896340735992468100)      /* 1/log(2)  */
10. #define EXP1       (2.71828182845904523536028747135)      /*   e       */
11. #define E          (2.71828182845904523536028747135)      /*   e       */
12. #define INV_E      (0.36787944117144232159552377016)      /*  e^-1     */
13. #define SQRT2      (1.41421356237309504880168872420)      /* sqrt(2)   */
14. #define INV_SQRT2  (0.70710678118654752440084436210)      /* sqrt(1/2) */
15. #define SQRT3      (1.73205080756887729352744634150)      /* sqrt(3)   */
16. #define INV_SQRT3  (0.57735026918962576450914878050)      /* sqrt(1/3) */
17.  
18.  
19. typedef struct
20. {
21.     double x;
22.     double y;
23. } complex;
24.  
25. double sin(double arg);
26. double cos(double arg);
27. double tan(double arg);
28.  
29. double asin(double arg);
30. double acos(double arg);
31. double atan(double arg);
32. double atan2(double arg1, double arg2);
33.  
34. double sinh(double arg);
35. double cosh(double arg);
36. double tanh(double arg);
37.  
38. double exp(double arg);
39. double log(double arg);
40. double log10(double arg);
41. double sqrt(double arg);
42.  
43. double floor(double x);
44. double ceil(double x);
45. double frexp(double x, int *exp);
46. double modf(double x, double *intptr);
47. double cabs(complex arg);
48.  
49. double hypot(double a, double b);
50. double gamma(double arg);
51. double erf(double arg);
52. double erfc(double arg);
53. double j0(double arg);
54. double y0(double arg);
55. double j1(double arg);
56. double y1(double arg);
57. double jn(int n, double x);
58. double yn(int n, double x);

Seno y coseno:
sin.c

Código: C

1. /*
2.   C program for floating point sin and cos.
3.   Calls modf.
4.   There are no error exits.
5.   Coefficients are #3370 from Hart & Cheney (18.80D).
6.  
7.   Error de los coeficientes (sin contar con el redondeo de las operaciones)
8.   error relativo = 1.5*sin(5*x) * 10^-19
9.   error absoluto < 1.5 * 10^-19
10.   62 bit de precision en la mantisa.
11. */
12. #include "cmath.h"
13.  
14. static const double p0 =  0.1357884097877375669092680e8;
15. static const double p1 = -0.4942908100902844161158627e7;
16. static const double p2 =  0.4401030535375266501944918e6;
17. static const double p3 = -0.1384727249982452873054457e5;
18. static const double p4 =  0.1459688406665768722226959e3;
19. static const double q0 =  0.8644558652922534429915149e7;
20. static const double q1 =  0.4081792252343299749395779e6;
21. static const double q2 =  0.9463096101538208180571257e4;
22. static const double q3 =  0.1326534908786136358911494e3;
23.  
24.  
25. static double sinus(double arg, int quad)
26. {
27.    double x, x2, temp;
28.    int integ;
29.  
30.    /* argumento negativo */
31.    if (arg < 0)
32.    {
33.       arg = -arg;
34.       quad = quad + 2;
35.    }
36.  
37.    /* argumento mayor de 90 grados */
38.    x = arg * TWO_PI;
39.    integ = (int) floor(x);
40.    x = x - integ;
41.  
42.    quad = (quad + integ) & 0x03;
43.    if (quad & 01)  x = 1.0 - x;
44.    if (quad > 1)   x =     - x;
45.  
46.    /* calcula el seno de x */
47.    x2 = x*x;
48.    temp = ((((p4*x2 + p3)*x2 + p2)*x2 + p1)*x2 + p0)*x;
49.    temp = temp/((((x2+q3)*x2+q2)*x2+q1)*x2+q0);
50.    return(temp);
51. }
52.  
53. double cos(double arg)
54. {
55.    if (arg<0)  arg = -arg;
56.    return(sinus(arg, 1));
57. }
58.  
59. double sin(double arg)
60. {
61.    return(sinus(arg, 0));
62. }

Arco seno y arco coseno
asin.c

Código: C

1. /*
2.    asin(arg) and acos(arg) return the arcsin, arccos,
3.    respectively of their arguments.
4.  
5.    Arctan is called after appropriate range reduction.
6. */
7. #include <errno.h>
8. #include "cmath.h"
9.  
10. double asin(double arg)
11. {
12.    double temp;
13.    int sign;
14.  
15.    sign = 1;
16.    if(arg < 0)
17.    {
18.       arg = -arg;
19.       sign = -1;
20.    }
21.  
22.    if(arg > 1.0)
23.    {
24.       errno = EDOM;
25.       return(0.0);
26.    }
27.  
28.    temp = sqrt(1.0 - arg*arg);
29.    if(arg > 0.7) temp = PI_TWO - atan(temp/arg);
30.    else          temp = atan(arg/temp);
31.  
32.    if (sign>0) return(temp);
33.    else        return(-temp);
34. }
35.  
36. double acos(double arg)
37. {
38.    if((arg > 1.0) || (arg < -1.0))
39.    {
40.       errno = EDOM;
41.       return(0.0);
42.    }
43.    return(PI_TWO - asin(arg));
44. }

Tangente
tan.c

Código: C

1. /*
2.    floating point tangent
3.    A series is used after range reduction.
4.    Coefficients are #4285 from Hart & Cheney. (19.74D)
5. */
6. #include <errno.h>
7. #include "cmath.h"
8.  
9. static const double p0 = -0.1306820264754825668269611177e+5;
10. static const double p1 =  0.1055970901714953193602353981e+4;
11. static const double p2 = -0.1550685653483266376941705728e+2;
12. static const double p3 =  0.3422554387241003435328470489e-1;
13. static const double p4 =  0.3386638642677172096076369000e-4;
14. static const double q0 = -0.1663895238947119001851464661e+5;
15. static const double q1 =  0.4765751362916483698926655581e+4;
16. static const double q2 = -0.1555033164031709966900124574e+3;
17.  
18. double tan(double arg)
19. {
20.    double temp, e, x, xx;
21.    int flag, sign;
22.  
23.    flag = 0;
24.    sign = 1;
25.    if(arg < 0.0)
26.    {
27.       arg = -arg;
28.       sign = -1;
29.    }
30.    arg = arg*INV_PI;   /*overflow?*/
31.    x = modf(arg, &e);
32.  
33.    switch((int)e & 0x03)
34.    {
35.       case 1: x = 1.0 - x;
36.               flag = 1;
37.               break;
38.  
39.       case 2: sign = - sign;
40.               flag = 1;
41.               break;
42.  
43.       case 3: x = 1.0 - x;
44.               sign = - sign;
45.               break;
46.  
47.       case 0: break;
48.    }
49.  
50.    xx = x*x;
51.    temp = ((((p4*xx+p3)*xx+p2)*xx+p1)*xx+p0)*x;
52.    temp = temp/(((xx+q2)*xx+q1)*xx+q0);
53.  
54.    if(flag == 1)
55.    {
56.       if(temp == 0.0)
57.       {
58.          errno = ERANGE;
59.          if (sign > 0) return(HUGE);
60.          else          return(-HUGE);
61.       }
62.       temp = 1.0/temp;
63.    }
64.    if (sign > 0) return(temp);
65.    else          return(-temp);
66. }

Arco tangente
atan.c

Código: C

1. /*
2.    floating-point arctangent
3.  
4.    atan returns the value of the arctangent of its
5.    argument in the range [-pi/2,pi/2].
6.  
7.    atan2 returns the arctangent of arg1/arg2
8.    in the range [-pi,pi].
9.  
10.    there are no error returns.
11.  
12.    coefficients are #5077 from Hart & Cheney. (19.56D)
13. */
14.  
15. static const double sq2p1 = 2.414213562373095048802e0;
16. static const double sq2m1 = 0.414213562373095048802e0;
17. static const double pio2  = 1.570796326794896619231e0;
18. static const double pio4  = 0.785398163397448309615e0;
19. static const double p4    = 0.161536412982230228262000000e2;
20. static const double p3    = 0.268425481955039737941410000e3;
21. static const double p2    = 0.115302935154048501154281360e4;
22. static const double p1    = 0.178040631643319697105464587e4;
23. static const double p0    = 0.896785974036638619599874880e3;
24. static const double q4    = 0.589569705084446222279100000e2;
25. static const double q3    = 0.536265374031215315104235000e3;
26. static const double q2    = 0.166678381488163371845217980e4;
27. static const double q1    = 0.207933497444540981287275926e4;
28. static const double q0    = 0.896785974036638619624811620e3;
29.  
30. static double xatan(double arg);
31. static double satan(double arg);
32.  
33.  
34. /*
35.    atan makes its argument positive and
36.    calls the inner routine satan.
37. */
38. double atan(double arg)
39. {
40.    
41.    if(arg>0) return(satan(arg));
42.    else      return(-satan(-arg));
43. }
44.  
45.  
46. /*
47.    atan2 discovers what quadrant the angle
48.    is in and calls atan.
49. */
50. double atan2(double arg1, double arg2)
51. {
52.    if((arg1+arg2)==arg1)
53.       if(arg1 >= 0.0) return(pio2);
54.       else return(-pio2);
55.    else if(arg2 < 0.0)
56.       if(arg1 >= 0.0) return(pio2+pio2 - satan(-arg1/arg2));
57.       else            return(-pio2-pio2 + satan(arg1/arg2));
58.    else if(arg1>0) return(satan(arg1/arg2));
59.    else            return(-satan(-arg1/arg2));
60. }
61.  
62. /*
63.    satan reduces its argument (known to be positive)
64.    to the range [0,0.414...] and calls xatan.
65. */
66.  
67. static double satan(double arg)
68. {
69.    if(arg < sq2m1)      return(xatan(arg));
70.    else if(arg > sq2p1) return(pio2 - xatan(1.0/arg));
71.    else                 return(pio4 + xatan((arg-1.0)/(arg+1.0)));
72. }
73.  
74. /*
75.    xatan evaluates a series valid in the
76.    range [-0.414...,+0.414...].
77. */
78. static double xatan(double arg)
79. {
80.    double argsq;
81.    double value;
82.  
83.    argsq = arg*arg;
84.    value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0);
85.    value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0);
86.    return(value*arg);
87. }

Exponencial:
exp.c

Código: C

1. /*
2.    exp returns the exponential function of its
3.    floating-point argument.
4.    The coefficients are #1069 from Hart and Cheney. (22.35D)
5. */
6. #include <errno.h>
7. #include "cmath.h"
8.  
9. extern int nerrno;
10. static const double  p0    = 0.2080384346694663001443843411e7;
11. static const double  p1    = 0.3028697169744036299076048876e5;
12. static const double  p2    = 0.6061485330061080841615584556e2;
13. static const double  q0    = 0.6002720360238832528230907598e7;
14. static const double  q1    = 0.3277251518082914423057964422e6;
15. static const double  q2    = 0.1749287689093076403844945335e4;
16. static const double  maxf  = 10000;
17.  
18. double exp(double arg)
19. {
20.    double x, xx, temp1, temp2;
21.    int ent;
22.  
23.    if(arg == 0.0)  return(1.0);
24.    if(arg < -maxf) return(0.0);
25.    if(arg > maxf)
26.    {
27.       errno = ERANGE;
28.       return(HUGE);
29.    }
30.  
31.    x = arg * INV_LN2;
32.    ent = (int)floor(x);
33.    x  = (x - ent) - 0.5;
34.    xx = x*x;
35.  
36.    temp1 = ((p2*xx+p1)*xx+p0)*x;
37.    temp2 = ((1.0*xx+q2)*xx+q1)*xx + q0;
38.    temp1 = SQRT2*(temp2+temp1)/(temp2-temp1);
39.    while(ent--) temp1 *= 2;
40.    return(temp1);
41. }

Logaritmo Natural
log.c

Código: C

1. /*
2.    log returns the natural logarithm of its floating
3.    point argument.
4.    The coefficients are #2705 from Hart & Cheney. (19.38D)
5.    It calls frexp.
6. */
7. #include <errno.h>
8. #include "cmath.h"
9.  
10. static const double p0 = -0.240139179559210510e2;
11. static const double p1 =  0.309572928215376501e2;
12. static const double p2 = -0.963769093368686593e1;
13. static const double p3 =  0.421087371217979714e0;
14. static const double q0 = -0.120069589779605255e2;
15. static const double q1 =  0.194809660700889731e2;
16. static const double q2 = -0.891110902798312337e1;
17.  
18. double log(double arg)
19. {
20.    double x,z, zz, temp;
21.    int exp;
22.  
23.    if(arg <= 0.0)
24.    {
25.       errno = EDOM;
26.       return(-HUGE);
27.    }
28.  
29.    x = frexp(arg, &exp);
30.  
31.    if(x < INV_SQRT2)
32.    {
33.       x *= 2;
34.       exp--;
35.    }
36.  
37.    z = (x-1)/(x+1);
38.    zz = z*z;
39.    temp = ((p3*zz + p2)*zz + p1)*zz + p0;
40.    temp = temp/(((1.0*zz + q2)*zz + q1)*zz + q0);
41.    temp = temp*z + exp*LN2;
42.    return(temp);
43. }
44.  
45. double log10(double arg)
46. {
47.    return(log(arg)*INV_LN10);
48. }

Potencias
pow.c

Código: C

1. /*
2.    computes a^b.
3.    uses log and exp
4. */
5. #include <errno.h>
6. #include "cmath.h"
7.  
8. double pow(double arg1, double arg2)
9. {
10.    double temp;
11.    long l;
12.  
13.    if(arg1 <= 0.0)
14.    {
15.       if(arg1 == 0.0 )
16.       {
17.          if(arg2 <= 0.0) errno = EDOM;
18.          return(0.0);
19.       }
20.       if(arg2 == 0.0) return(1.0);
21.  
22.       l = (long)arg2;
23.       if((double)l != arg2)
24.       {
25.          errno = EDOM;
26.          return(0.0);
27.       }
28.  
29.       temp = exp(arg2 * log(-arg1));
30.       if(l & 1L) return(-temp);
31.       return(temp);
32.    }
33.    return(exp(arg2 * log(arg1)));
34. }

Raiz cuadrada:
sqrt.c

Código: C

1. /*
2.         sqrt returns the square root of its floating
3.         point argument. Newton's method.
4.  
5.         calls frexp
6. */
7. #include <errno.h>
8. #include "cmath.h"
9.  
10. double sqrt(double arg)
11. {
12.    double x, temp;
13.    int exp, i;
14.  
15.    if(arg <= 0.0)
16.    {
17.       if(arg < 0.0) errno = EDOM;
18.       return(0.0);
19.    }
20.  
21.    x = frexp(arg, &exp);
22.    if(exp & 1)
23.    {
24.       x *= 2;
25.       exp--;
26.    }
27.    
28.    temp = 0.5*(1.0 + x);
29.  
30.    while(exp > 60)
31.    {
32.       temp *= (1L<<30);
33.       exp -= 60;
34.    }
35.    if(exp >= 0) temp *= 1L << (exp/2);
36.  
37.    while(exp < -60)
38.    {
39.       temp /= (1L<<30);
40.       exp += 60;
41.    }
42.    if(exp < 0)  temp /= 1L << (-exp/2);
43.  
44.    i = 5;
45.    do {
46.       if (i-- == 0) break;
47.       x = temp;
48.       temp = 0.5*(temp + arg/temp);
49.    } while (temp != x);
50.  
51.    return(temp);
52. }

Seno y coseno hiperbólico
sinh.c

Código: C

1. /*
2.    sinh(arg) returns the hyperbolic sine of its floating-
3.    point argument.
4.  
5.    The exponential function is called for arguments
6.    greater in magnitude than 0.5.
7.  
8.    A series is used for arguments smaller in magnitude than 0.5.
9.    The coefficients are #2029 from Hart & Cheney. (20.36D)
10.  
11.    cosh(arg) is computed from the exponential function for
12.    all arguments.
13. */
14. #include "cmath.h"
15.  
16. static const double p0  = -0.6307673640497716991184787251e+6;
17. static const double p1  = -0.8991272022039509355398013511e+5;
18. static const double p2  = -0.2894211355989563807284660366e+4;
19. static const double p3  = -0.2630563213397497062819489e+2;
20. static const double q0  = -0.6307673640497716991212077277e+6;
21. static const double q1   = 0.1521517378790019070696485176e+5;
22. static const double q2  = -0.173678953558233699533450911e+3;
23.  
24. double sinh(double arg)
25. {
26.    double temp, argsq;
27.    int sign;
28.  
29.    sign = 1;
30.    if(arg < 0)
31.    {
32.       arg = - arg;
33.       sign = -1;
34.    }
35.  
36.    if(arg > 21.0)
37.    {
38.       temp = exp(arg)/2;
39.       if (sign>0) return(temp);
40.       else        return(-temp);
41.    }
42.    
43.    if(arg > 0.5)
44.    {
45.       if (sign>0) return((exp(arg) - exp(-arg))/2);
46.       else        return((-exp(arg) - exp(-arg))/2);
47.    }
48.    
49.    argsq = arg*arg;
50.    temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg;
51.    temp /= (((argsq+q2)*argsq+q1)*argsq+q0);
52.  
53.    if (sign>0) return(temp);
54.    else        return(-temp);
55. }
56.  
57. double cosh(double arg)
58. {
59.    if(arg < 0)     arg = - arg;
60.    if(arg > 21.0)  return(exp(arg)/2);
61.  
62.    return((exp(arg) + exp(-arg))/2);
63. }

Tangente hiperbólica
tanh.c

Código: C

1. /*
2.    tanh(arg) computes the hyperbolic tangent of its floating
3.    point argument.
4.  
5.    sinh and cosh are called except for large arguments, which
6.    would cause overflow improperly.
7. */
8.  
9. double tanh(double arg)
10. {
11.    int sign;
12.  
13.    sign = 1;
14.    if(arg < 0.)
15.    {
16.       arg = -arg;
17.       sign = -1;
18.    }
19.  
20.    if(arg > 21.0)
21.    {
22.       if (sign>0) return(1.0);
23.       else        return(-1.0);
24.    }
25.  
26.    if (sign>0) return(sinh(arg)/cosh(arg));
27.    else        return(sinh(arg)/cosh(arg));
28. }

Saludos.

[Picuino]

Biblioteca en C de las funciones matemáticas más comunes en precisión doble (64 bits):

Continuación con funciones de Bessel:

Según Wikipedia

Aplicaciones
La Ecuación de Bessel aparece cuando se buscan soluciones a la ecuación de Laplace o a la ecuación de Helmholtz por el método de separación de variables en coordenadas cilíndricas o esféricas. Por ello, las funciones de Bessel son especialmente importantes en muchos problemas de propagación de ondas, potenciales estáticos y cualquier otro problema descrito por las ecuaciones de Helmholtz o Laplace en simetrías cilíndricas o esféricas.

Función de Bessel de orden 0:
j0.c

Código: C

1. /*
2.    floating point Bessel's function
3.    of the first and second kinds
4.    of order zero
5.  
6.    j0(x) returns the value of J0(x)
7.    for all real values of x.
8.  
9.    There are no error returns.
10.  
11.    There is a niggling bug in J0 which
12.    causes errors up to 2e-16 for x in the
13.    interval [-8,8].
14.    The bug is caused by an inappropriate order
15.    of summation of the series.  rhm will fix it
16.    someday.
17.  
18.    Coefficients are from Hart & Cheney.
19.    #5849 (19.22D)
20.    #6549 (19.25D)
21.    #6949 (19.41D)
22.  
23.    y0(x) returns the value of Y0(x)
24.    for positive real values of x.
25.    For x<=0, error number EDOM is set and a
26.    large negative value is returned.
27.  
28.    The values of Y0 have not been checked
29.    to more than ten places.
30.  
31.    Coefficients are from Hart & Cheney.
32.    #6245 (18.78D)
33.    #6549 (19.25D)
34.    #6949 (19.41D)
35. */
36.  
37. #include <errno.h>
38. #include "cmath.h"
39.  
40. static void asympt(double arg);
41.  
42. static double pzero, qzero;
43. static const double tpi  = 0.6366197723675813430755350535e0;
44. static const double pio4 = 0.7853981633974483096156608458e0;
45. static const double p1[] = {
46.     0.4933787251794133561816813446e21,
47.    -0.1179157629107610536038440800e21,
48.     0.6382059341072356562289432465e19,
49.    -0.1367620353088171386865416609e18,
50.     0.1434354939140344111664316553e16,
51.    -0.8085222034853793871199468171e13,
52.     0.2507158285536881945555156435e11,
53.    -0.4050412371833132706360663322e8,
54.     0.2685786856980014981415848441e5,
55. };
56. static const double q1[] = {
57.    0.4933787251794133562113278438e21,
58.    0.5428918384092285160200195092e19,
59.    0.3024635616709462698627330784e17,
60.    0.1127756739679798507056031594e15,
61.    0.3123043114941213172572469442e12,
62.    0.66999876729822396718140286600e9,
63.    0.11146360984629853781824025430e7,
64.    0.13630636523289706044428105070e4,
65.    1.0
66. };
67. static const double p2[] = {
68.    0.5393485083869438325262122897e7,
69.    0.1233238476817638145232406055e8,
70.    0.8413041456550439208464315611e7,
71.    0.2016135283049983642487182349e7,
72.    0.1539826532623911470917825993e6,
73.    0.2485271928957404011288128951e4,
74.    0.0,
75. };
76. static const double q2[] = {
77.    0.5393485083869438325560444960e7,
78.    0.1233831022786324960844856182e8,
79.    0.8426449050629797331554404810e7,
80.    0.2025066801570134013891035236e7,
81.    0.1560017276940030940592769933e6,
82.    0.2615700736920839685159081813e4,
83.    1.0,
84. };
85. static const double p3[] = {
86.    -0.3984617357595222463506790588e4,
87.    -0.1038141698748464093880530341e5,
88.    -0.8239066313485606568803548860e4,
89.    -0.2365956170779108192723612816e4,
90.    -0.2262630641933704113967255053e3,
91.    -0.4887199395841261531199129300e1,
92.     0.0,
93. };
94.  
95. static const double q3[] = {
96.    0.2550155108860942382983170882e6,
97.    0.6667454239319826986004038103e6,
98.    0.5332913634216897168722255057e6,
99.    0.1560213206679291652539287109e6,
100.    0.1570489191515395519392882766e5,
101.    0.4087714673983499223402830260e3,
102.    1.0,
103. };
104. static const double p4[] = {
105.    -.2750286678629109583701933175e20,
106.    0.6587473275719554925999402049e20,
107.    -.5247065581112764941297350814e19,
108.    0.1375624316399344078571335453e18,
109.    -.1648605817185729473122082537e16,
110.    0.1025520859686394284509167421e14,
111.    -.3436371222979040378171030138e11,
112.    0.5915213465686889654273830069e8,
113.    -.4137035497933148554125235152e5,
114. };
115. static const double q4[] = {
116.    0.3726458838986165881989980e21,
117.    0.4192417043410839973904769661e19,
118.    0.2392883043499781857439356652e17,
119.    0.9162038034075185262489147968e14,
120.    0.2613065755041081249568482092e12,
121.    0.5795122640700729537480087915e9,
122.    0.1001702641288906265666651753e7,
123.    0.1282452772478993804176329391e4,
124.    1.0,
125. };
126.  
127. double j0(double arg)
128. {
129.    double argsq, n, d;
130.    int i;
131.  
132.    if(arg < 0.0) arg = -arg;
133.    if(arg > 8.0)
134.    {
135.       asympt(arg);
136.       n = arg - pio4;
137.       return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
138.    }
139.    argsq = arg*arg;
140.  
141.    for(n=0,d=0,i=8;i>=0;i--)
142.    {
143.       n = n*argsq + p1[i];
144.       d = d*argsq + q1[i];
145.    }
146.    return(n/d);
147. }
148.  
149. double y0(double arg)
150. {
151.    double argsq, n, d;
152.    int i;
153.  
154.    errno = 0;
155.    if(arg <= 0.0)
156.    {
157.       errno = EDOM;
158.       return(-HUGE);
159.    }
160.    if(arg > 8.)
161.    {
162.       asympt(arg);
163.       n = arg - pio4;
164.       return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
165.    }
166.    argsq = arg*arg;
167.    for(n=0,d=0,i=8;i>=0;i--)
168.    {
169.       n = n*argsq + p4[i];
170.       d = d*argsq + q4[i];
171.    }
172.    return(n/d + tpi*j0(arg)*log(arg));
173. }
174.  
175. static void asympt(double arg)
176. {
177.    double zsq, n, d;
178.    int i;
179.    zsq = 64.0/(arg*arg);
180.    for(n=0,d=0,i=6;i>=0;i--)
181.    {
182.       n = n*zsq + p2[i];
183.       d = d*zsq + q2[i];
184.    }
185.    pzero = n/d;
186.    for(n=0,d=0,i=6;i>=0;i--)
187.    {
188.       n = n*zsq + p3[i];
189.       d = d*zsq + q3[i];
190.    }
191.    qzero = (8.0/arg)*(n/d);
192. }

Función de Bessel de orden 1:
j1.c

Código: C

1. /*
2.    floating point Bessel's function
3.    of the first and second kinds
4.    of order one
5.  
6.    j1(x) returns the value of J1(x)
7.    for all real values of x.
8.  
9.    There are no error returns.
10.  
11.    There is a niggling bug in J1 which
12.    causes errors up to 2e-16 for x in the
13.    interval [-8,8].
14.    The bug is caused by an inappropriate order
15.    of summation of the series.  rhm will fix it
16.    someday.
17.  
18.    Coefficients are from Hart & Cheney.
19.    #6050 (20.98D)
20.    #6750 (19.19D)
21.    #7150 (19.35D)
22.  
23.    y1(x) returns the value of Y1(x)
24.    for positive real values of x.
25.    For x<=0, error number EDOM is set and a
26.    large negative value is returned.
27.  
28.    The values of Y1 have not been checked
29.    to more than ten places.
30.  
31.    Coefficients are from Hart & Cheney.
32.    #6447 (22.18D)
33.    #6750 (19.19D)
34.    #7150 (19.35D)
35. */
36. #include <errno.h>
37. #include "cmath.h"
38.  
39. static void asympt(double arg);
40.  
41. static double pzero, qzero;
42. static const double tpi = .6366197723675813430755350535e0;
43. static const double pio4   = .7853981633974483096156608458e0;
44. static const double p1[] = {
45.    0.581199354001606143928050809e21,
46.    -.6672106568924916298020941484e20,
47.    0.2316433580634002297931815435e19,
48.    -.3588817569910106050743641413e17,
49.    0.2908795263834775409737601689e15,
50.    -.1322983480332126453125473247e13,
51.    0.3413234182301700539091292655e10,
52.    -.4695753530642995859767162166e7,
53.    0.2701122710892323414856790990e4,
54. };
55. static const double q1[] = {
56.    0.1162398708003212287858529400e22,
57.    0.1185770712190320999837113348e20,
58.    0.6092061398917521746105196863e17,
59.    0.2081661221307607351240184229e15,
60.    0.5243710262167649715406728642e12,
61.    0.1013863514358673989967045588e10,
62.    0.1501793594998585505921097578e7,
63.    0.1606931573481487801970916749e4,
64.    1.0,
65. };
66. static const double p2[] = {
67.    -.4435757816794127857114720794e7,
68.    -.9942246505077641195658377899e7,
69.    -.6603373248364939109255245434e7,
70.    -.1523529351181137383255105722e7,
71.    -.1098240554345934672737413139e6,
72.    -.1611616644324610116477412898e4,
73.    0.0,
74. };
75. static const double q2[] = {
76.    -.4435757816794127856828016962e7,
77.    -.9934124389934585658967556309e7,
78.    -.6585339479723087072826915069e7,
79.    -.1511809506634160881644546358e7,
80.    -.1072638599110382011903063867e6,
81.    -.1455009440190496182453565068e4,
82.    1.0,
83. };
84. static const double p3[] = {
85.    0.3322091340985722351859704442e5,
86.    0.8514516067533570196555001171e5,
87.    0.6617883658127083517939992166e5,
88.    0.1849426287322386679652009819e5,
89.    0.1706375429020768002061283546e4,
90.    0.3526513384663603218592175580e2,
91.    0.0,
92. };
93. static const double q3[] = {
94.    0.7087128194102874357377502472e6,
95.    0.1819458042243997298924553839e7,
96.    0.1419460669603720892855755253e7,
97.    0.4002944358226697511708610813e6,
98.    0.3789022974577220264142952256e5,
99.    0.8638367769604990967475517183e3,
100.    1.0,
101. };
102. static const double p4[] = {
103.    -.9963753424306922225996744354e23,
104.    0.2655473831434854326894248968e23,
105.    -.1212297555414509577913561535e22,
106.    0.2193107339917797592111427556e20,
107.    -.1965887462722140658820322248e18,
108.    0.9569930239921683481121552788e15,
109.    -.2580681702194450950541426399e13,
110.    0.3639488548124002058278999428e10,
111.    -.2108847540133123652824139923e7,
112.    0.0,
113. };
114. static const double q4[] = {
115.    0.5082067366941243245314424152e24,
116.    0.5435310377188854170800653097e22,
117.    0.2954987935897148674290758119e20,
118.    0.1082258259408819552553850180e18,
119.    0.2976632125647276729292742282e15,
120.    0.6465340881265275571961681500e12,
121.    0.1128686837169442121732366891e10,
122.    0.1563282754899580604737366452e7,
123.    0.1612361029677000859332072312e4,
124.    1.0,
125. };
126.  
127. double j1(double arg)
128. {
129.    double xsq, n, d, x;
130.    int i;
131.  
132.    x = arg;
133.    if(x < 0.0) x = -x;
134.    if(x > 8.0)
135.    {
136.       asympt(x);
137.       n = x - 3.*pio4;
138.       n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
139.       if(arg <0.) n = -n;
140.       return(n);
141.    }
142.    xsq = x*x;
143.    for(n=0,d=0,i=8;i>=0;i--)
144.    {
145.       n = n*xsq + p1[i];
146.       d = d*xsq + q1[i];
147.    }
148.    return(arg*n/d);
149. }
150.  
151. double y1(double arg)
152. {
153.    double xsq, n, d, x;
154.    int i;
155.  
156.    errno = 0;
157.    x = arg;
158.    if(x <= 0.0)
159.    {
160.       errno = EDOM;
161.       return(-HUGE);
162.    }
163.    if(x > 8.0)
164.    {
165.       asympt(x);
166.       n = x - 3*pio4;
167.       return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
168.    }
169.    xsq = x*x;
170.    for(n=0,d=0,i=9;i>=0;i--)
171.    {
172.       n = n*xsq + p4[i];
173.       d = d*xsq + q4[i];
174.    }
175.    return(x*n/d + tpi*(j1(x)*log(x)-1.0/x));
176. }
177.  
178. static void asympt(double arg)
179. {
180.    double zsq, n, d;
181.    int i;
182.    zsq = 64./(arg*arg);
183.    for(n=0,d=0,i=6;i>=0;i--)
184.    {
185.       n = n*zsq + p2[i];
186.       d = d*zsq + q2[i];
187.    }
188.    pzero = n/d;
189.    for(n=0,d=0,i=6;i>=0;i--)
190.    {
191.       n = n*zsq + p3[i];
192.       d = d*zsq + q3[i];
193.    }
194.    qzero = (8.0/arg)*(n/d);
195. }

Función de Bessel de orden n:
jn.c

Código: C

1. /*
2.    floating point Bessel's function of
3.    the first and second kinds and of
4.    integer order.
5.  
6.    int n;
7.    double x;
8.    jn(n,x);
9.  
10.    returns the value of Jn(x) for all
11.    integer values of n and all real values
12.    of x.
13.  
14.    There are no error returns.
15.    Calls j0, j1.
16.  
17.    For n=0, j0(x) is called,
18.    for n=1, j1(x) is called,
19.    for n<x, forward recursion us used starting
20.    from values of j0(x) and j1(x).
21.    for n>x, a continued fraction approximation to
22.    j(n,x)/j(n-1,x) is evaluated and then backward
23.    recursion is used starting from a supposed value
24.    for j(n,x). The resulting value of j(0,x) is
25.    compared with the actual value to correct the
26.    supposed value of j(n,x).
27.  
28.    yn(n,x) is similar in all respects, except
29.    that forward recursion is used for all
30.    values of n>1.
31. */
32. #include <errno.h>
33. #include "cmath.h"
34.  
35.  
36. double jn(int n, double x)
37. {
38.    int i;
39.    double a, b, temp;
40.    double xsq, t;
41.  
42.    if(n<0)
43.    {
44.       n = -n;
45.       x = -x;
46.    }
47.    if(n==0) return(j0(x));
48.    if(n==1) return(j1(x));
49.    if(x == 0.0) return(0.);
50.  
51.    if(n<=(int)x)
52.    {
53.       a = j0(x);
54.       b = j1(x);
55.       for(i=1;i<n;i++)
56.       {
57.          temp = b;
58.          b = (2.0*i/x)*b - a;
59.          a = temp;
60.       }
61.       return(b);
62.    }
63.  
64.    xsq = x*x;
65.    for(t=0,i=n+16;i>n;i--) t = xsq/(2.0*i - t);
66.    t = x/(2.0*n-t);
67.    a = t;
68.    b = 1;
69.    for(i=n-1; i>0; i--)
70.    {
71.       temp = b;
72.       b = (2.0*i/x)*b - a;
73.       a = temp;
74.    }
75.    return(t*j0(x)/b);
76. }
77.  
78. double yn(int n, double x)
79. {
80.    int i, sign;
81.    double a, b, temp;
82.  
83.    if (x <= 0)
84.    {
85.       errno = EDOM;
86.       return(-HUGE);
87.    }
88.    sign = 1;
89.    if(n<0)
90.    {
91.       n = -n;
92.       if(n%2 == 1) sign = -1;
93.    }
94.    if(n==0) return(y0(x));
95.    if(n==1)
96.    {
97.       if (sign<0) return(-y1(x));
98.       return(y1(x));
99.    }
100.  
101.    a = y0(x);
102.    b = y1(x);
103.    for(i=1;i<n;i++)
104.    {
105.       temp = b;
106.       b = (2.0*i/x)*b - a;
107.       a = temp;
108.    }
109.  
110.    if (sign<0) return(-b);
111.    return(b);
112. }

Saludos.

[Picuino]

Biblioteca en C de funciones matemáticas en precisión doble (64 bits):

Función Gamma: Calcula el factorial de números enteros o reales.


Función error:  Calcula la distribución Normal.




Función Error
erf.c

Código: C

1. /*
2.    C program for floating point error function
3.  
4.    erf(x) returns the error function of its argument
5.    erfc(x) returns 1.0-erf(x)
6.  
7.    erf(x) is defined by
8.    ${2 over sqrt(pi)} int from 0 to x e sup {-t sup 2} dt$
9.  
10.    the entry for erfc is provided because of the
11.    extreme loss of relative accuracy if erf(x) is
12.    called for large x and the result subtracted
13.    from 1. (e.g. for x= 10, 12 places are lost).
14.  
15.    There are no error returns.
16.    Coefficients for large x are #5667 from Hart & Cheney (18.72D).
17. */
18. #include <errno.h>
19. #include "cmath.h"
20.  
21. #define M 7
22. #define N 9
23.  
24. static const double torp = 1.1283791670955125738961589031;
25. static const double p1[] = {
26.     0.804373630960840172832162e5,
27.     0.740407142710151470082064e4,
28.     0.301782788536507577809226e4,
29.     0.380140318123903008244444e2,
30.     0.143383842191748205576712e2,
31.    -0.288805137207594084924010e0,
32.     0.007547728033418631287834e0,
33. };
34. static const double q1[]  = {
35.     0.804373630960840172826266e5,
36.     0.342165257924628539769006e5,
37.     0.637960017324428279487120e4,
38.     0.658070155459240506326937e3,
39.     0.380190713951939403753468e2,
40.     0.100000000000000000000000e1,
41.     0.0,
42. };
43. static const double p2[]  = {
44.     0.1826334884229511259216899900e4,
45.     0.2898029329216765561127584600e4,
46.     0.2320439590251635247384768711e4,
47.     0.1143262070703886173606073338e4,
48.     0.3685196154710010637133875746e3,
49.     0.7708161730368428609781633646e2,
50.     0.9675807882987265400604202961e1,
51.     0.5641877825507397413087057563e0,
52.     0.0,
53. };
54. static const double q2[]  = {
55.     0.1826334884229511259557643800e4,
56.     0.4958827564721140714954384220e4,
57.     0.6089542423272443550463306800e4,
58.     0.4429612803883682726711528526e4,
59.     0.2094384367789539593790281779e4,
60.     0.6617361207107653469211984771e3,
61.     0.1371255960500622202878443578e3,
62.     0.1714980943627607849376131193e2,
63.     1.0,
64. };
65.  
66. double erf(double arg)
67. {
68.    int sign, i;
69.    double argsq, d, n;
70.  
71.    sign = 1;
72.    if(arg < 0.0)
73.    {
74.       arg = -arg;
75.       sign = -1;
76.    }
77.    if(arg < 0.5)
78.    {
79.       argsq = arg*arg;
80.       for(n=0,d=0,i=M-1; i>=0; i--)
81.       {
82.          n = n*argsq + p1[i];
83.          d = d*argsq + q1[i];
84.       }
85.       if (sign<0) return(-torp*arg*n/d);
86.       return(torp*arg*n/d);
87.    }
88.    if(arg >= 10.0)  return(sign*1.);
89.    if (sign<0) return(-1.0 + erfc(arg));
90.    return(1.0 - erfc(arg));
91. }
92.  
93. double erfc(double arg)
94. {
95.    double n, d;
96.    int i;
97.  
98.    errno = 0;
99.    if(arg < 0.0) return(2.0 - erfc(-arg));
100.    if(arg < 0.5) return(1.0 - erf(arg));
101.    if(arg >= 10.0) return(0.0);
102.  
103.    for(n=0,d=0,i=N-1; i>=0; i--)
104.    {
105.       n = n*arg + p2[i];
106.       d = d*arg + q2[i];
107.    }
108.    return(exp(-arg*arg)*n/d);
109. }

Función Gamma:
gamma.c

Código: C

1. /*
2.    C program for floating point log gamma function
3.  
4.    gamma(x) computes the log of the absolute
5.    value of the gamma function.
6.    The sign of the gamma function is returned in the
7.    external quantity signgam.
8.  
9.    The coefficients for expansion around zero
10.    are #5243 from Hart & Cheney; for expansion
11.    around infinity they are #5404.
12. */
13. #include <errno.h>
14. #include "cmath.h"
15.  
16. int   signgam = 0;
17. static const double goobie = 0.9189385332046727417803297;
18. static const double pi     = 3.1415926535897932384626434;
19. #define M 6
20. #define N 8
21.  
22. static const double p1[] = {
23.     0.83333333333333101837e-1,
24.    -0.277777777735865004e-2,
25.     0.793650576493454e-3,
26.    -0.5951896861197e-3,
27.     0.83645878922e-3,
28.    -0.1633436431e-2,
29. };
30. static const double p2[] = {
31.    -0.42353689509744089647e5,
32.    -0.20886861789269887364e5,
33.    -0.87627102978521489560e4,
34.    -0.20085274013072791214e4,
35.    -0.43933044406002567613e3,
36.    -0.50108693752970953015e2,
37.    -0.67449507245925289918e1,
38.     0.0,
39. };
40. static const double q2[] = {
41.    -0.42353689509744090010e5,
42.    -0.29803853309256649932e4,
43.     0.99403074150827709015e4,
44.    -0.15286072737795220248e4,
45.    -0.49902852662143904834e3,
46.     0.18949823415702801641e3,
47.    -0.23081551524580124562e2,
48.     0.10000000000000000000e1,
49. };
50.  
51. static double pos(double arg);
52. static double neg(double arg);
53. static double asym(double arg);
54.  
55.  
56. double gamma(double arg)
57. {
58.    signgam = 1.0;
59.    if (arg <= 0.0) return(neg(arg));
60.    if (arg > 8.0)  return(asym(arg));
61.    return(log(pos(arg)));
62. }
63.  
64. static double asym(double arg)
65. {
66.    double n, argsq;
67.    int i;
68.  
69.    argsq = 1./(arg*arg);
70.    for(n=0,i=M-1; i>=0; i--)
71.       n = n*argsq + p1[i];
72.  
73.    return((arg-0.5)*log(arg) - arg + goobie + n/arg);
74. }
75.  
76. static double neg(double arg)
77. {
78.    double temp;
79.  
80.    arg = -arg;
81.    temp = sin(pi*arg);
82.    if(temp == 0.0)
83.    {
84.       errno = EDOM;
85.       return(HUGE);
86.    }
87.    if(temp < 0.) temp = -temp;
88.    else signgam = -1;
89.    return(-log(arg*pos(arg)*temp/pi));
90. }
91.  
92. static double pos(double arg)
93. {
94.    double n, d, s;
95.    register i;
96.  
97.    if(arg < 2.0) return(pos(arg+1.0)/arg);
98.    if(arg > 3.0) return((arg-1.0)*pos(arg-1.0));
99.  
100.    s = arg - 2.0;
101.    for(n=0,d=0,i=N-1; i>=0; i--)
102.    {
103.       n = n*s + p2[i];
104.       d = d*s + q2[i];
105.    }
106.    return(n/d);
107. }

Saludos.