Analizaremos dos aproximaciones polinómicas: La serie de Taylor y la aproximación de Padé (esta última se trata realmente de una aproximación racional).
Serie de Taylor
Lo primero que se le suele venir a uno a la cabeza cuando piensa en aproximaciones polinómicas suele ser la serie de Taylor, dicha serie es muy sencilla de calcular y genera una buena aproximación en el entorno de un punto. En este caso se ha optado por aproximar la función exponencial en el entorno de x=0:
$$e^{x} \simeq 1 + {x \over 1!} + {x^2 \over 2!} + {x^3 \over 3!} + ...$$
Esto nos da una muy buena aproximación, aunque para conseguir un error aceptable, es necesario calcular la serie de Taylor para un orden relativamente alto. El error con respecto a la función “exp” de la librería matemática de C comienza a ser asumible a partir del orden 6 trabajando en coma flotante.
Aproximación de Padé
La aproximación de Padé de orden (m, n) es la función racional:
$$R(x) = {p_0 + p_{1}x + p_{2}x^2 + ... + p_{m}x^m \over 1 + q_{1}x+ q_{2}x^2 + ... + q_{n}x^n}$$
Que cumple que:
$$f(0) = R(0)$$
$$f'(0) = R'(0)$$
$$f''(0) = R''(0)$$
$$...$$
$$f^{(m+n)}(0) = R^{(m+n)}(0)$$
El cálculo de los coeficientes de Padé no es trivial y existen varias técnicas para obtenerlos, como el algoritmo Epsilon de Wynn o el algoritmo euclídeo extendido para el cálculo del máximo común divisor. Por suerte, para la función exponencial, podemos consultar las tablas de Padé que ya se encuentran en internet calculadas para diferentes órdenes (valores de m y de n):
Aquí una tabla publicada en wikipedia.
Se ha optado en este caso por utilizar la aproximación de Padé de orden [3 / 3] (m = n = 3).
A continuación puede verse una implementación en C de ambas aproximaciones.
double taylor_exp(double x) { double ret = 0; int i; double num = 1; double den = 1; for (i = 0; i <= 6; i++) { ret += num / den; num *= x; den *= (i + 1); } return ret; } double pade_exp(double x) { double x2 = x * x; double x3 = x2 * x; double num = 1 + (x / 2) + (x2 / 10) + (x3 / 120); double den = 1 - (x / 2) + (x2 / 10) - (x3 / 120); return num / den; }
Como puede apreciarse, la serie de Taylor es de orden 6, mientras que la aproximación de Padé que se ha implementado es la de orden [3 / 3]. A continuación se reproduce la salida de una prueba de ambas funciones comparándolas con la función “exp” de la librería matemática:
# ./taylor_vs_pade_float
exp(0.250000):
exp() function : 1.2840254167
6th order taylor : 1.2840254042
3rd order pade : 1.2840254175
Como puede verse, la aproximación de Padé consigue un error comparable al de la serie de Taylor con muchas menos operaciones.
Utilizar aritmética de punto fijo
Ahora que están ambos algoritmos implementados en coma flotante, pasaremos el cálculo a aritmética de punto fijo en formato Q16.16 (más info sobre la notación Q). En el formato Q16.16 tenemos 16 bits para la parte entera y 16 bits para la parte fraccionaria, en total 32 bits.
typedef int32_t fixedpoint_t; #define FP_NEG(x) (-(x)) #define FP_ADD(x, y) ((x) + (y)) #define FP_SUB(x, y) ((x) - (y)) #define FP_MUL(x, y) ((int32_t) (((int64_t) (x)) * ((int64_t) (y)) >> 16)) #define FP_DIV(x, y) ((int32_t) ((((int64_t) (x)) << 16) / ((int64_t) (y)))) #define TO_FP(x) ((int32_t) ((x) << 16)) #define FROM_FP(x) ((x) >> 16) #define FP_FRACTIONAL_BITS 16
Las funciones anteriores puede ser ahora reescritas para utilizar el formato Q16.16:
fixedpoint_t taylor_exp(fixedpoint_t x) { fixedpoint_t ret = 0; int i; fixedpoint_t num = TO_FP(1); fixedpoint_t den = TO_FP(1); for (i = 0; i <= 6; i++) { ret = FP_ADD(ret, FP_DIV(num, den)); num = FP_MUL(num, x); den = FP_MUL(den, FP_ADD(TO_FP(i), TO_FP(1))); } return ret; } fixedpoint_t pade_exp(fixedpoint_t x) { fixedpoint_t x2 = FP_MUL(x, x); fixedpoint_t x3 = FP_MUL(x2, x); fixedpoint_t num = FP_ADD(TO_FP(1), FP_ADD(FP_DIV(x, TO_FP(2)), FP_ADD(FP_DIV(x2, TO_FP(10)), FP_DIV(x3, TO_FP(120))))); fixedpoint_t den = FP_ADD(FP_SUB(FP_DIV(x2, TO_FP(10)), FP_DIV(x3, TO_FP(120))), FP_SUB(TO_FP(1), FP_DIV(x, TO_FP(2)))); return FP_DIV(num, den); }
En este caso, realizando la misma prueba obtenemos resultados algo peores (debido a la pérdida de precisión inherente al uso del punto fijo) y, aunque para las dos aproximaciones obtenemos el mismo valor, la aproximación de Padé requiere menor cantidad de operaciones que la serie de Taylor.
# ./taylor_vs_pade_fixed
exp(0.250000):
exp() function : 1.2840254167
6th order taylor : 1.2839965820
3rd order pade : 1.2839965820
La aproximación de Padé, como ha podido verse, da mejores resultados que las series de Taylor como aproximación a la función exponencial, tanto desde el punto de vista de la eficiencia como de la precisión.
[ añadir comentario ] ( 3417 visualizaciones ) | [ 0 trackbacks ] | enlace permanente |
( 3 / 2166 )
Calendario
