
/* 10^x.  Hastings again */

#include <math.h>

#define N_COEFF 7

static double _expt_coeff[7] =
	{
	1.15129277603,
	0.66273088429,
	0.25439357484,
	0.07295173666,
	0.01742111988,
	0.00255491796,
	0.00093264267
	};

/* this stuff should probly get split out into someplace else... */
#define N_PWRS_OF_TEN 16
int _max_pwr_of_ten = N_PWRS_OF_TEN;
double _pwrs_of_ten[N_PWRS_OF_TEN + 1] =
	{
	1.0,
	10.0,
	100.0,
	1000.0,
	10000.0,
	100000.0,
	1000000.0,
	10000000.0,
	100000000.0,
	1000000000.0,
	10000000000.0,
	100000000000.0,
	1000000000000.0,
	10000000000000.0,
	100000000000000.0,
	1000000000000000.0,
	10000000000000000.0,
	};

double expt(x)
double x;
{
  int offset, invert, i;
  double accum, x0, xn;
  
  x0 = x;
  if (x0 == 0.0)
	return(1.0);
  if (x0 < 0.0)
  	{
	invert = 1;
	x0 = -x0;
	}
    else
	invert = 0;
  offset = x0;			/* get int part */
  x0 = x0 - offset;		/* leave fractional part */
  xn = x0;
  accum = 1.0;
  for (i = 0 ; i < N_COEFF ; i++)
	{
	accum = accum + (xn * _expt_coeff[i]);
	xn = xn * x0;
	}
  accum = accum * accum;	/* almost forgot that last 2 */
  for ( ; offset > 0 ; )
	{
	i = (offset > _max_pwr_of_ten) ? _max_pwr_of_ten : offset;
	accum = accum * _pwrs_of_ten[i];
	offset -= i;
	}
  if (invert)
	accum = 1.0 / accum;
  return(accum);
}


#ifdef TEST
main()
{
  double x;
  
  for (x = -2.0 ; x < 3.0 ; x += .2)
	printf("%f: %f\n", x, expt(x));
}
#endif
