pow_zi.c   [plain text]

#include "f2c.h"

extern void z_div (doublecomplex *, doublecomplex *, doublecomplex *);
void
pow_zi (doublecomplex * p, doublecomplex * a, integer * b)	/* p = a**b  */
{
  integer n;
  unsigned long u;
  double t;
  doublecomplex q, x;
  static doublecomplex one = { 1.0, 0.0 };

  n = *b;
  q.r = 1;
  q.i = 0;

  if (n == 0)
    goto done;
  if (n < 0)
    {
      n = -n;
      z_div (&x, &one, a);
    }
  else
    {
      x.r = a->r;
      x.i = a->i;
    }

  for (u = n;;)
    {
      if (u & 01)
	{
	  t = q.r * x.r - q.i * x.i;
	  q.i = q.r * x.i + q.i * x.r;
	  q.r = t;
	}
      if (u >>= 1)
	{
	  t = x.r * x.r - x.i * x.i;
	  x.i = 2 * x.r * x.i;
	  x.r = t;
	}
      else
	break;
    }
done:
  p->i = q.i;
  p->r = q.r;
}