1 /* 2 * "@(#)pow_zi.c 1.1" 3 */ 4 5 #include "complex" 6 7 pow_zi(p, a, b) /* p = a**b */ 8 dcomplex *p, *a; 9 long int *b; 10 { 11 long int n; 12 double t; 13 dcomplex x; 14 15 n = *b; 16 p->dreal = 1; 17 p->dimag = 0; 18 19 if(n == 0) 20 return; 21 if(n < 0) 22 { 23 n = -n; 24 z_div(&x, p, a); 25 } 26 else 27 { 28 x.dreal = a->dreal; 29 x.dimag = a->dimag; 30 } 31 32 for( ; ; ) 33 { 34 if(n & 01) 35 { 36 t = p->dreal * x.dreal - p->dimag * x.dimag; 37 p->dimag = p->dreal * x.dimag + p->dimag * x.dreal; 38 p->dreal = t; 39 } 40 if(n >>= 1) 41 { 42 t = x.dreal * x.dreal - x.dimag * x.dimag; 43 x.dimag = 2 * x.dreal * x.dimag; 44 x.dreal = t; 45 } 46 else 47 break; 48 } 49 } 50