/* powi.c - calc x^n, where n is an integer! */

/* Very slightly modified version of power() from Computer Language, Sept. 86,
	pg 87, by Jon Snader (who extracted the binary algorithm from Donald Knuth,
	"The Art of Computer Programming", vol 2, 1969).
	powi() will only be called when an exponentiation with an integer
	exponent is performed, thus tests & code for fractional exponents were 
	removed.	-- cw lightfoot, COBALT BLUE
  COBALT BLUE has NO copyright on the functions in this file.  The source may 
  be considered public domain.
 */
#include "powi.h"
#define IS_ODD(j)	((j) & 1 )

/* this will only define powi if this is not an alpha running the native os */
#if !defined(__alpha) || defined(__linux)

/* #ifndef __alpha   this fcn is on an alpha, only defn if not alpha */

double powi( double x, long int n )	/* returns:  x^n */
/* x;	 base */
/* n;	 exponent */
{
	double p;	/* holds partial product */

	if( x == 0 )
		return( 0. );

	/* test for negative exponent */
	if( n < 0 )
	{	
		n = -n;
		x = 1/x;
	}

	p = IS_ODD(n) ? x : 1;	/* test & set zero power */

	while( n >>= 1 ){	/* now do the other powers */
		x *= x;			/* sq previous power of x */
		if( IS_ODD(n) )	/* if low order bit set */
			p *= x;		/*	then, multiply partial product by latest power of x */
		}

	return( p );
}

#	endif /*alpha unix*/

	/*		*		*		*		*/

long ipow( long m, long n )	/* returns:  m^n */
/* m;		 base */
/* n;		 exponent */
{
long p;	/* holds partial product */

if( m == 0 || (n < 0 && m > 1) )
	return( 0L );
	/* NOTE: negative exponent always results in 0 for integers!
	 * (except for the case when m==1 or -1) */

if( n < 0 ){	/* test for negative exponent */
	n = -n;
	m = 1/m;
	}

p = IS_ODD(n) ? m : 1;	/* test & set zero power */

while( n >>= 1 ){	/* now do the other powers */
	m *= m;			/* sq previous power of m */
	if( IS_ODD(n) )	/* if low order bit set */
		p *= m;		/*	then, multiply partial product by latest power of m */
	}

return( p );
}