1#include "f2c.h" 2 3#ifdef KR_headers 4VOID pow_zi(p, a, b) /* p = a**b */ 5 doublecomplex *p, *a; integer *b; 6#else 7extern void z_div(doublecomplex*, doublecomplex*, doublecomplex*); 8void pow_zi(doublecomplex *p, doublecomplex *a, integer *b) /* p = a**b */ 9#endif 10{ 11 integer n; 12 unsigned long u; 13 double t; 14 doublecomplex q, x; 15 static doublecomplex one = {1.0, 0.0}; 16 17 n = *b; 18 q.r = 1; 19 q.i = 0; 20 21 if(n == 0) 22 goto done; 23 if(n < 0) 24 { 25 n = -n; 26 z_div(&x, &one, a); 27 } 28 else 29 { 30 x.r = a->r; 31 x.i = a->i; 32 } 33 34 for(u = n; ; ) 35 { 36 if(u & 01) 37 { 38 t = q.r * x.r - q.i * x.i; 39 q.i = q.r * x.i + q.i * x.r; 40 q.r = t; 41 } 42 if(u >>= 1) 43 { 44 t = x.r * x.r - x.i * x.i; 45 x.i = 2 * x.r * x.i; 46 x.r = t; 47 } 48 else 49 break; 50 } 51 done: 52 p->i = q.i; 53 p->r = q.r; 54 } 55