110516Sdlw /* 222939Skre * Copyright (c) 1980 Regents of the University of California. 322939Skre * All rights reserved. The Berkeley software License Agreement 422939Skre * specifies the terms and conditions for redistribution. 522939Skre * 6*33391Sbostic * @(#)pow_zi.c 5.2 01/24/88 710516Sdlw */ 810516Sdlw 910516Sdlw #include "complex" 1010516Sdlw 11*33391Sbostic #define Z_MULEQ(A,B) \ 12*33391Sbostic t = (A).dreal * (B).dreal - (A).dimag * (B).dimag,\ 13*33391Sbostic (A).dimag = (A).dreal * (B).dimag + (A).dimag * (B).dreal,\ 14*33391Sbostic (A).dreal = t /* A *= B */ 15*33391Sbostic 16*33391Sbostic void 1710516Sdlw pow_zi(p, a, b) /* p = a**b */ 18*33391Sbostic dcomplex *p, *a; 19*33391Sbostic long int *b; 2010516Sdlw { 21*33391Sbostic register long n = *b; 22*33391Sbostic double t; 23*33391Sbostic dcomplex x; 2410516Sdlw 25*33391Sbostic x = *a; 26*33391Sbostic p->dreal = (double)1, p->dimag = (double)0; 27*33391Sbostic if (!n) 28*33391Sbostic return; 29*33391Sbostic if (n < 0) { 30*33391Sbostic z_div(&x, p, a); 31*33391Sbostic n = -n; 3210516Sdlw } 33*33391Sbostic while (!(n&1)) { 34*33391Sbostic Z_MULEQ(x, x); 35*33391Sbostic n >>= 1; 3610516Sdlw } 37*33391Sbostic for (*p = x; --n > 0; Z_MULEQ(*p, x)) 38*33391Sbostic while (!(n&1)) { 39*33391Sbostic Z_MULEQ(x, x); 40*33391Sbostic n >>= 1; 4110516Sdlw } 4210516Sdlw } 43