Aproximación de Padé para el cálculo eficiente de la función exponencial 
Implementar la función exponencial en un sistema embebido con poca RAM, poca memoria de programa y sin coprocesador matemático pasa, normalmente, por intentar evitar el uso de la librería matemática de C. La sobrecarga que produciría el utilizar la función “exp” de dicha librería unida a la sobrecarga propia de la manipulación de datos en coma flotante por software desaconsejan totalmente el uso de dicha librería en sistema embebidos pequeños. Analizaremos diferentes aproximaciones polinomiales a la función exponencial y el uso de aritmética de punto fijo para realizar dicho cálculo.

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.

Comentarios 
Lo sentimos. No se permiten nuevos comentarios después de 90 días.