| #include "f2c.h" |
| |
| #ifdef KR_headers |
| VOID pow_zi(p, a, b) /* p = a**b */ |
| doublecomplex *p, *a; integer *b; |
| #else |
| extern void z_div(doublecomplex*, doublecomplex*, doublecomplex*); |
| void pow_zi(doublecomplex *p, doublecomplex *a, integer *b) /* p = a**b */ |
| #endif |
| { |
| 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; |
| } |