51 lines
737 B
C
51 lines
737 B
C
#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;
|
|
}
|