/*Copyright (C) 1992, 1996 by Thomas Glen Smith.  All Rights Reserved.*/
/* power APL2 V1.0.0 ***************************************************
* Called by powerx.                                                    *
* Power is a front-end to the standard library function pow().  X and Y*
* are its two arguments, and it returns the Yth power of X, according  *
* to the following rules:                                              *
* 1.  if absolute value of Y is less than FUZZ, return 1.0 regardless  *
*     of the value of x.                                               *
* 2.  X may be less than zero if Y is an approximation of a rational   *
*     number with an odd denominator.                                  *
***********************************************************************/
#define INCLUDES MATH
#include "includes.h"
double power(left,rite)
double left,rite;
{
	Dabs; Mod; Pow;
	extern int aplerr;
	extern double fuzz;
	double absdif,absleft,absrite,cur,den,dif,num,ret;
	int i;

	absrite = (rite >= 0.0) ? rite : -rite; /* absolute value */
	if (absrite <= fuzz) return(1.0);
	absleft = (left >= 0.0) ? left : -left; /* absolute value */
	if (mod(absrite,1.0) > fuzz)
		ret=pow(absleft,rite);
	else { /* right is approximate integer */
		i=absrite; /* convert to int. */
		ret=absleft;
		while(0<--i) ret *= absleft;
		if (rite < 0.0) ret = 1.0 / ret;
	}
	if (left < 0.0) { /* determine if denominator is odd */
		num = den = 1.0; /* start with 1/1 */
		for ( i = 0; i < 400; i++ ) {
			dif = rite - (cur = num / den);
			absdif = (dif >= 0.0) ? dif : -dif; /* absolute value */
			if (absdif <= fuzz) {
				if (mod(den,2.0) < fuzz) aplerr = 86;
				if (mod(num,2.0) < fuzz) return(ret);
				return(-ret);
			}
			if (rite > cur) num += 1.0;
			else den += 1.0;
		}
		aplerr = 87; /* bad exponent for negative number */
	}
	return(ret);
}