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