/* Copyright (C) 1992 by Thomas Glen Smith.  All Rights Reserved. */
/* powerx APL2 V1.0.0  *************************************************
* Power for complex numbers.                                           *
***********************************************************************/
#define INCLUDES MATH+TRIGKEYS
#include "includes.h"
#define ABS(v) ((v) > 0e0 ? (v) : -(v))
void powerx(left,rite,ret)
double *left,*rite,*ret;
{
	Circular; Circulax; Dabsx; Dividex; Mod;
	Power; Powerx; Timesx;
     extern int aplerr;
	extern double fuzz;
	double absval,n,neg,wrka[2],wrkb[2],wrkc[2];
     int i;

	
	if (ABS(*(rite+1)) > fuzz) { /* rite complex? */
     	*wrka = *rite;
          *(wrka+1) = 0e0;
          powerx(left,wrka,wrkb); /* wrkb = left * a */
          *wrka = circular(COS,*(left+1));
          *(wrka+1) = circular(SIN,*(left+1)); /*wrka=cos(b)+iXsin(b)*/
          timesx(wrkb,wrka,ret);
          return;
     }
	/* control reaches here only if rite is a real number */
	if (*(left+1) <= fuzz && *(left+1) >= fuzz) { /* left real? */
     	*(ret+1) = 0e0;
          *ret = power(*left,*rite);
          return;
     }
     else { /* left is complex */
	     neg = (*rite < 0e0) ? -1.0 : 1.0;
          absval = neg * *rite;
		if (mod(absval, 1.0) > fuzz) { /* Is it non-integer power? */
			/* a*1%n = magnitude:(|a}*1%n, angle:alpha%n */
			n = 1.0 / absval;
			if (mod(n, 1.0) > fuzz) { /* non-integer n in 1%n? */
				aplerr = 996; /* non-integer powers of complex numbers */
		          return;
			}
			dabsx(left,wrka); /* Get magnitude: |a. */
			*wrkb = power(*wrka,absval); /* Get (|a)*1%n. */
			*(wrkb+1) = atan2(*(left+1),*left) / n;
			*wrka = *wrkb * cos(*(wrkb+1));
			*(wrka + 1) = *wrkb * sin(*(wrkb+1));
	     } else {
			/* right is approximate integer */
	          i = absval; /* absolute integer value of rite */
	          XFR(wrka,left) /* copy left to wrka */
	          while(--i) {
	          	timesx(wrka,left,wrkb);
	          	XFR(wrka,wrkb) /* copy wrkb to wrka */
	          }
		}
          if (neg == -1.0) {
          	*wrkb = 1.0;
               *(wrkb+1) = 0.0;
          	dividex(wrkb,wrka,ret);
		}
          else XFR(ret,wrka);
     }
}