WebSVN - shark - Blame - Rev 2 - /shark/branches/pj/oslib/libm/msun/src/e

Rev	Author	Line No.	Line
2	pj	1	/* @(#)e_exp.c 5.1 93/09/24 */
		2	/*
		3	* ====================================================
		4	* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
		5	*
		6	* Developed at SunPro, a Sun Microsystems, Inc. business.
		7	* Permission to use, copy, modify, and distribute this
		8	* software is freely granted, provided that this notice
		9	* is preserved.
		10	* ====================================================
		11	*/
		12
		13	#ifndef lint
		14	static char rcsid[] = "$\Id: e_exp.c,v 1.3.2.1 1997/02/23 11:03:02 joerg Exp $";
		15	#endif
		16
		17	/* __ieee754_exp(x)
		18	* Returns the exponential of x.
		19	*
		20	* Method
		21	* 1. Argument reduction:
		22	* Reduce x to an r so that \|r\| <= 0.5*ln2 ~ 0.34658.
		23	* Given x, find r and integer k such that
		24	*
		25	* x = kln2 + r, \|r\| <= 0.5ln2.
		26	*
		27	* Here r will be represented as r = hi-lo for better
		28	* accuracy.
		29	*
		30	* 2. Approximation of exp(r) by a special rational function on
		31	* the interval [0,0.34658]:
		32	* Write
		33	* R(r*2) = r(exp(r)+1)/(exp(r)-1) = 2 + rr/6 - r*4/360 + ...
		34	* We use a special Reme algorithm on [0,0.34658] to generate
		35	* a polynomial of degree 5 to approximate R. The maximum error
		36	* of this polynomial approximation is bounded by 2**-59. In
		37	* other words,
		38	* R(z) ~ 2.0 + P1z + P2z*2 + P3z*3 + P4z*4 + P5z**5
		39	* (where z=r*r, and the values of P1 to P5 are listed below)
		40	* and
		41	* \| 5 \| -59
		42	* \| 2.0+P1z+...+P5z - R(z) \| <= 2
		43	* \| \|
		44	* The computation of exp(r) thus becomes
		45	* 2*r
		46	* exp(r) = 1 + -------
		47	* R - r
		48	* r*R1(r)
		49	* = 1 + r + ----------- (for better accuracy)
		50	* 2 - R1(r)
		51	* where
		52	* 2 4 10
		53	* R1(r) = r - (P1r + P2r + ... + P5*r ).
		54	*
		55	* 3. Scale back to obtain exp(x):
		56	* From step 1, we have
		57	* exp(x) = 2^k * exp(r)
		58	*
		59	* Special cases:
		60	* exp(INF) is INF, exp(NaN) is NaN;
		61	* exp(-INF) is 0, and
		62	* for finite argument, only exp(0)=1 is exact.
		63	*
		64	* Accuracy:
		65	* according to an error analysis, the error is always less than
		66	* 1 ulp (unit in the last place).
		67	*
		68	* Misc. info.
		69	* For IEEE double
		70	* if x > 7.09782712893383973096e+02 then exp(x) overflow
		71	* if x < -7.45133219101941108420e+02 then exp(x) underflow
		72	*
		73	* Constants:
		74	* The hexadecimal values are the intended ones for the following
		75	* constants. The decimal values may be used, provided that the
		76	* compiler will convert from decimal to binary accurately enough
		77	* to produce the hexadecimal values shown.
		78	*/
		79
		80	#include "math.h"
		81	#include "math_private.h"
		82
		83	#ifdef __STDC__
		84	static const double
		85	#else
		86	static double
		87	#endif
		88	one = 1.0,
		89	halF[2] = {0.5,-0.5,},
		90	huge = 1.0e+300,
		91	twom1000= 9.33263618503218878990e-302, /* 2*-1000=0x01700000,0/
		92	o_threshold= 7.09782712893383973096e+02, /* 0x40862E42, 0xFEFA39EF */
		93	u_threshold= -7.45133219101941108420e+02, /* 0xc0874910, 0xD52D3051 */
		94	ln2HI[2] ={ 6.93147180369123816490e-01, /* 0x3fe62e42, 0xfee00000 */
		95	-6.93147180369123816490e-01,},/* 0xbfe62e42, 0xfee00000 */
		96	ln2LO[2] ={ 1.90821492927058770002e-10, /* 0x3dea39ef, 0x35793c76 */
		97	-1.90821492927058770002e-10,},/* 0xbdea39ef, 0x35793c76 */
		98	invln2 = 1.44269504088896338700e+00, /* 0x3ff71547, 0x652b82fe */
		99	P1 = 1.66666666666666019037e-01, /* 0x3FC55555, 0x5555553E */
		100	P2 = -2.77777777770155933842e-03, /* 0xBF66C16C, 0x16BEBD93 */
		101	P3 = 6.61375632143793436117e-05, /* 0x3F11566A, 0xAF25DE2C */
		102	P4 = -1.65339022054652515390e-06, /* 0xBEBBBD41, 0xC5D26BF1 */
		103	P5 = 4.13813679705723846039e-08; /* 0x3E663769, 0x72BEA4D0 */
		104
		105
		106	#ifdef __STDC__
		107	double __generic___ieee754_exp(double x) /* default IEEE double exp */
		108	#else
		109	double __generic___ieee754_exp(x) /* default IEEE double exp */
		110	double x;
		111	#endif
		112	{
		113	double y,hi=0.0,lo=0.0,c,t;
		114	int32_t k=0,xsb;
		115	u_int32_t hx;
		116
		117	GET_HIGH_WORD(hx,x);
		118	xsb = (hx>>31)&1; /* sign bit of x */
		119	hx &= 0x7fffffff; /* high word of \|x\| */
		120
		121	/* filter out non-finite argument */
		122	if(hx >= 0x40862E42) { /* if \|x\|>=709.78... */
		123	if(hx>=0x7ff00000) {
		124	u_int32_t lx;
		125	GET_LOW_WORD(lx,x);
		126	if(((hx&0xfffff)\|lx)!=0)
		127	return x+x; /* NaN */
		128	else return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
		129	}
		130	if(x > o_threshold) return hugehuge; / overflow */
		131	if(x < u_threshold) return twom1000twom1000; / underflow */
		132	}
		133
		134	/* argument reduction */
		135	if(hx > 0x3fd62e42) { /* if \|x\| > 0.5 ln2 */
		136	if(hx < 0x3FF0A2B2) { /* and \|x\| < 1.5 ln2 */
		137	hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
		138	} else {
		139	k = invln2*x+halF[xsb];
		140	t = k;
		141	hi = x - tln2HI[0]; / tln2HI is exact here /
		142	lo = t*ln2LO[0];
		143	}
		144	x = hi - lo;
		145	}
		146	else if(hx < 0x3e300000) { /* when \|x\|<2*-28 /
		147	if(huge+x>one) return one+x;/* trigger inexact */
		148	}
		149	else k = 0;
		150
		151	/* x is now in primary range */
		152	t = x*x;
		153	c = x - t(P1+t(P2+t(P3+t(P4+t*P5))));
		154	if(k==0) return one-((x*c)/(c-2.0)-x);
		155	else y = one-((lo-(x*c)/(2.0-c))-hi);
		156	if(k >= -1021) {
		157	u_int32_t hy;
		158	GET_HIGH_WORD(hy,y);
		159	SET_HIGH_WORD(y,hy+(k<<20)); /* add k to y's exponent */
		160	return y;
		161	} else {
		162	u_int32_t hy;
		163	GET_HIGH_WORD(hy,y);
		164	SET_HIGH_WORD(y,hy+((k+1000)<<20)); /* add k to y's exponent */
		165	return y*twom1000;
		166	}
		167	}

Subversion Repositories shark

(root)/shark/branches/pj/oslib/libm/msun/src/e_exp.c - Rev 2