110516Sdlw /* 2*22939Skre * Copyright (c) 1980 Regents of the University of California. 3*22939Skre * All rights reserved. The Berkeley software License Agreement 4*22939Skre * specifies the terms and conditions for redistribution. 5*22939Skre * 6*22939Skre * @(#)pow_zi.c 5.1 06/07/85 710516Sdlw */ 810516Sdlw 910516Sdlw #include "complex" 1010516Sdlw 1110516Sdlw pow_zi(p, a, b) /* p = a**b */ 1210516Sdlw dcomplex *p, *a; 1310516Sdlw long int *b; 1410516Sdlw { 1510516Sdlw long int n; 1610516Sdlw double t; 1710516Sdlw dcomplex x; 1810516Sdlw 1910516Sdlw n = *b; 2010516Sdlw p->dreal = 1; 2110516Sdlw p->dimag = 0; 2210516Sdlw 2310516Sdlw if(n == 0) 2410516Sdlw return; 2510516Sdlw if(n < 0) 2610516Sdlw { 2710516Sdlw n = -n; 2810516Sdlw z_div(&x, p, a); 2910516Sdlw } 3010516Sdlw else 3110516Sdlw { 3210516Sdlw x.dreal = a->dreal; 3310516Sdlw x.dimag = a->dimag; 3410516Sdlw } 3510516Sdlw 3610516Sdlw for( ; ; ) 3710516Sdlw { 3810516Sdlw if(n & 01) 3910516Sdlw { 4010516Sdlw t = p->dreal * x.dreal - p->dimag * x.dimag; 4110516Sdlw p->dimag = p->dreal * x.dimag + p->dimag * x.dreal; 4210516Sdlw p->dreal = t; 4310516Sdlw } 4410516Sdlw if(n >>= 1) 4510516Sdlw { 4610516Sdlw t = x.dreal * x.dreal - x.dimag * x.dimag; 4710516Sdlw x.dimag = 2 * x.dreal * x.dimag; 4810516Sdlw x.dreal = t; 4910516Sdlw } 5010516Sdlw else 5110516Sdlw break; 5210516Sdlw } 5310516Sdlw } 54