dpowe.c

#include "mathutils.h"

#define LN2 (.6931471806)
#define C (68243.0)
#define TOL .00001
#define RUNAWAY 100


/*****************************************************************
	A *quick and dirty* function for approximating the natural exponent
	(e^y) for a double floating-point precision real y.
	Based entirely off of Schraudolph's findings, 1999. C likewise
	tuned based off of '99 paper.
  ***************************************************************/
double fast_dpowe(double y){
	double y_a = ABS(y);	
	split_d y_split;
	y_split.bits.j = 0;
	
	double a = ((1.0*(0x1 << 20))/LN2);
	double b = (0x3FF << 20)*1.0;
	y_split.bits.i = (a*y_a) + (b - C);
	
	if(y < 0.0)
		return (1/y_split.d);
	return y_split.d;
}
/*****************************************************************
	Classic Taylor Series approximation for computing e^y for a 
	real precision point value y. Attempts to halt on convergence 
	for a specified tolerance and stops after 20 iterations to 
	avoid divergence
  ***************************************************************/
double taylor_dpowe(double y){
	double y_poly,e_y,fact,r,next,abs_y;
	abs_y = ABS(y); 
	y_poly = fact = 1;
	e_y = 0.0;
	int _i = 1;
	do{
		y_poly *= abs_y;
		fact *= _i;
		next = y_poly/fact;
		e_y += next;
		r = next/e_y;
		_i++;
	}while(r > TOL && _i < RUNAWAY);
	e_y = (y < 0) ? (1.0/(e_y + 1.0)) : 1.0 + e_y;
	return e_y;
}