/* * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands. * See the copyright notice in the ACK home directory, in the file "Copyright". * * Author: Ceriel J.H. Jacobs *//* $Header: /cvsup/minix/src/lib/ack/libp/exp.c,v 1.1 2005/10/10 15:27:46 beng Exp $ */#define __NO_DEFS#include <math.h>#include <pc_err.h>extern	_trp();#if __STDC__#include <float.h>#include <pc_math.h>#define M_MIN_D	DBL_MIN#define M_MAX_D	DBL_MAX#define M_DMINEXP DBL_MIN_EXP#endif#undef HUGE#define HUGE	1e1000static doubleLdexp(fl,exp)	double fl;	int exp;{	extern double _fef();	int sign = 1;	int currexp;	if (fl<0) {		fl = -fl;		sign = -1;	}	fl = _fef(fl,&currexp);	exp += currexp;	if (exp > 0) {		while (exp>30) {			fl *= (double) (1L << 30);			exp -= 30;		}		fl *= (double) (1L << exp);	}	else	{		while (exp<-30) {			fl /= (double) (1L << 30);			exp += 30;		}		fl /= (double) (1L << -exp);	}	return sign * fl;}double_exp(x)	double x;{	/*	Algorithm and coefficients from:			"Software manual for the elementary functions"			by W.J. Cody and W. Waite, Prentice-Hall, 1980	*/	static double p[] = {	        0.25000000000000000000e+0,	        0.75753180159422776666e-2,		0.31555192765684646356e-4	};	static double q[] = {	        0.50000000000000000000e+0,	        0.56817302698551221787e-1,	        0.63121894374398503557e-3,		0.75104028399870046114e-6	};	double	xn, g;	int	n;	int	negative = x < 0;	if (x <= M_LN_MIN_D) {		g = M_MIN_D/4.0;		if (g != 0.0) {			/* unnormalized numbers apparently exist */			if (x < (M_LN2 * (M_DMINEXP - 53))) return 0.0;		}		else {			if (x < M_LN_MIN_D) return 0.0;			return M_MIN_D;		}	}	if (x >= M_LN_MAX_D) {		if (x > M_LN_MAX_D) {			_trp(EEXP);			return HUGE;		}		return M_MAX_D;	}	if (negative) x = -x;	n = x * M_LOG2E + 0.5;	/* 1/ln(2) = log2(e), 0.5 added for rounding */	xn = n;	{		double	x1 = (long) x;		double	x2 = x - x1;		g = ((x1-xn*0.693359375)+x2) - xn*(-2.1219444005469058277e-4);	}	if (negative) {		g = -g;		n = -n;	}	xn = g * g;	x = g * POLYNOM2(xn, p);	n += 1;	return (Ldexp(0.5 + x/(POLYNOM3(xn, q) - x), n));}