/*
exp returns exponential of floating point argument

the coefficients are #1069 from Hart and Cheney. (22.35D)
*/

static	double	p0	= 0.2080384346694663001443843411e7;
static	double	p1	= 0.3028697169744036299076048876e5;
static	double	p2	= 0.6061485330061080804165584556e2;
static	double	q0	= 0.6002720360238832528230907598e7;
static	double	q1	= 0.3277251518082914423057964422e6;
static	double	q2	= 0.1749287689093076403844945335e4;
static	double	log2e	= 1.4426950408889634073599247;
static	double	sqrt2	= 1.4142135623730950488016887;

double	exp(arg)

	double	arg;
{
	double	fract;
	double	temp1,temp2,xsq;
	int	ent;
extern	double	floor();
extern	double	ldexp();
extern	double	itof();

	if (arg == 0.0) return(1.0);
	if (arg < -10000.0) return(0.0);
	if (arg > 10000.0) return(.99999999e35);

	arg = arg*log2e;
	ent = ftoi(floor(arg));
	fract = (arg-itof(ent))-0.5;
	xsq = fract*fract;
	temp1 = ((p2*xsq+p1)*xsq+p0)*fract;
	temp2 = ((xsq+q2)*xsq+q1)*xsq+q0;
	return(ldexp(sqrt2*(temp2+temp1)/(temp2-temp1),ent));
}