/*Copyright (C) 1992, 1996 by Thomas Glen Smith.  All Rights Reserved.*/
/* powerp APL2 V1.0.0  *************************************************
* 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"
#define ABS(v) ((v) > 0e0 ? (v) : -(v))
void powerp(left,rite,pret)
double *left,*rite,*pret; /* Warning - result may be a complex number. */
{
	Mod; Pow; Powerpx;
	extern int aplerr;
	extern double fuzz;
	double absdif,absleft,absrite,cur,den,dif,num,ret[2];
	int i;

	*(ret + 1) = 0e0; /* Default result is not complex. */
	absrite = ABS(*rite); /* absolute value */
	if (absrite <= fuzz) { *pret =(1.0); return; }
	absleft = ABS(*left); /* absolute value */
	if (mod(absrite,1.0) > fuzz) /* Is rite other than integer? */
		*ret = pow(absleft,*rite); /* Yes. */
	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 when left neg.*/
		num = den = 1.0; /* start with 1/1 */
		for ( i = 0; i < 400; i++ ) { /* num / den is est. of exponent. */
			cur = num / den; /* Estimate of exponent. */
			absdif = ABS(*rite - cur); /* Error of estimate. */
			if (absdif <= fuzz) { /* Estimate close enough? */
				if (mod(den,2.0) < fuzz) /* Denominator even? */
					powerpx(left,rite,pret,ret,num,den); /*Complex*/
				else if (mod(num,2.0) < fuzz) /* Numerator even? */
					*pret = *ret;
				else *pret = -*ret;
				*(pret+1) = *(ret+1);
                    return;
			}
			if (*rite > cur)
				num += 1.0; /* Revise estimate bigger. */
			else den += 1.0; /* Revise Estimate smaller. */
		}
		aplerr = 87; /* bad exponent for negative number */
	}
	*pret = *ret;
	*(pret+1) = *(ret+1);
}