text: /*
   * $Id: math_expf.c,v 1.3 2006-01-08 12:04:23 clib2devs Exp $
   *
  
   *
   *
   * PowerPC math library based in part on work by Sun Microsystems
   * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   *
   * Developed at SunPro, a Sun Microsystems, Inc. business.
   * Permission to use, copy, modify, and distribute this
   * software is freely granted, provided that this notice
   * is preserved.
   *
   *
   * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
   */
  
  #ifndef _MATH_HEADERS_H
  #include "math_headers.h"
  #endif /* _MATH_HEADERS_H */
  
  static const float
  	one = 1.0,
  	halF[2] = {
  		0.5,
  		-0.5,
  },
  	huge = 1.0e+30, twom100 = 7.8886090522e-31, /* 2**-100=0x0d800000 */
  	o_threshold = 8.872283911e+01,				/* 0x42b17218 */
  	u_threshold = -1.0397208405e+02,			/* 0xc2cff1b5 */
  	ln2HI[2] = {
  		6.9313812256e-01, /* 0x3f317180 */
  		-6.9313812256e-01,
  }, /* 0xbf317180 */
  	ln2LO[2] = {
  		9.0580006145e-06, /* 0x3717f7d1 */
  		-9.0580006145e-06,
  },							   /* 0xb717f7d1 */
  	invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
  	P1 = 1.6666667163e-01,	   /* 0x3e2aaaab */
  	P2 = -2.7777778450e-03,	   /* 0xbb360b61 */
  	P3 = 6.6137559770e-05,	   /* 0x388ab355 */
  	P4 = -1.6533901999e-06,	   /* 0xb5ddea0e */
  	P5 = 4.1381369442e-08;	   /* 0x3331bb4c */
  
  float expf(float x)
  {
  	float y, hi = 0.0, lo = 0.0, c, t;
  	LONG k = 0, xsb;
  	ULONG hx;
  
  	GET_FLOAT_WORD(hx, x);
  	xsb = (hx >> 31) & 1; /* sign bit of x */
  	hx &= 0x7fffffff;	  /* high word of |x| */
  
  	/* filter out non-finite argument */
  	if (hx >= 0x42b17218)
  	{ /* if |x|>=88.721... */
  		if (hx == 0x42b17218)
  		{ /* if |x| == 88.7228... */
  			if (xsb == 0)
  				return FLT_MAX;
  		}
  		if (hx > 0x7f800000)
  			return x + x; /* NaN */
  		if (hx == 0x7f800000)
  			return (xsb == 0) ? x : 0.0; /* exp(+-inf)={inf,0} */
  		if (x > o_threshold)
  			return huge * huge; /* overflow */
  		if (x < u_threshold)
  			return twom100 * twom100; /* underflow */
  	}
  
  	/* argument reduction */
  	if (hx > 0x3eb17218)
  	{ /* if  |x| > 0.5 ln2 */
  		if (hx < 0x3F851592)
  		{ /* and |x| < 1.5 ln2 */
  			hi = x - ln2HI[xsb];
  			lo = ln2LO[xsb];
  			k = 1 - xsb - xsb;
  		}
  		else
  		{
  			k = invln2 * x + halF[xsb];
  			t = k;
  			hi = x - t * ln2HI[0]; /* t*ln2HI is exact here */
  			lo = t * ln2LO[0];
  		}
  		x = hi - lo;
  	}
  	else if (hx < 0x31800000)
  	{ /* when |x|<2**-28 */
  		if (huge + x > one)
  			return one + x; /* trigger inexact */
  	}
  	else
  		k = 0;
  
  	/* x is now in primary range */
  	t = x * x;
  	c = x - t * (P1 + t * (P2 + t * (P3 + t * (P4 + t * P5))));
  	if (k == 0)
  		return one - ((x * c) / (c - (float)2.0) - x);
  	else
  		y = one - ((lo - (x * c) / ((float)2.0 - c)) - hi);
  	if (k >= -125)
  	{
  		ULONG hy;
  		GET_FLOAT_WORD(hy, y);
  		SET_FLOAT_WORD(y, hy + (k << 23)); /* add k to y's exponent */
  		return y;
  	}
  	else
  	{
  		ULONG hy;
  		GET_FLOAT_WORD(hy, y);
  		SET_FLOAT_WORD(y, hy + ((k + 100) << 23)); /* add k to y's exponent */
  		return y * twom100;
  	}
  }
  