110510Sdlw /* 222926Skre * Copyright (c) 1980 Regents of the University of California. 322926Skre * All rights reserved. The Berkeley software License Agreement 422926Skre * specifies the terms and conditions for redistribution. 522926Skre * 6*33391Sbostic * @(#)pow_ci.c 5.3 01/24/88 710510Sdlw */ 810510Sdlw 910510Sdlw #include "complex" 1010510Sdlw 11*33391Sbostic #ifdef tahoe 12*33391Sbostic 13*33391Sbostic #define C_MULEQ(A,B) \ 14*33391Sbostic t = (A).real * (B).real - (A).imag * (B).imag,\ 15*33391Sbostic (A).imag = (A).real * (B).imag + (A).imag * (B).real,\ 16*33391Sbostic (A).real = t /* A *= B */ 17*33391Sbostic 18*33391Sbostic void 1910510Sdlw pow_ci(p, a, b) /* p = a**b */ 20*33391Sbostic complex *p, *a; 21*33391Sbostic long *b; 2210510Sdlw { 23*33391Sbostic register long n = *b; 24*33391Sbostic register float t; 25*33391Sbostic complex x; 2610510Sdlw 27*33391Sbostic x = *a; 28*33391Sbostic p->real = (float)1, p->imag = (float)0; 29*33391Sbostic if (!n) 30*33391Sbostic return; 31*33391Sbostic if (n < 0) { 32*33391Sbostic c_div(&x, p, a); 33*33391Sbostic n = -n; 34*33391Sbostic } 35*33391Sbostic while (!(n&1)) { 36*33391Sbostic C_MULEQ(x, x); 37*33391Sbostic n >>= 1; 38*33391Sbostic } 39*33391Sbostic for (*p = x; --n > 0; C_MULEQ(*p, x)) 40*33391Sbostic while (!(n&1)) { 41*33391Sbostic C_MULEQ(x, x); 42*33391Sbostic n >>= 1; 43*33391Sbostic } 44*33391Sbostic } 4510510Sdlw 46*33391Sbostic #else /* !tahoe */ 4710510Sdlw 48*33391Sbostic extern void pow_zi(); 4929965Smckusick 50*33391Sbostic void 5129965Smckusick pow_ci(p, a, b) /* p = a**b */ 52*33391Sbostic complex *p, *a; 53*33391Sbostic long *b; 5429965Smckusick { 55*33391Sbostic dcomplex p1, a1; 5629965Smckusick 57*33391Sbostic a1.dreal = a->real; 58*33391Sbostic a1.dimag = a->imag; 5929965Smckusick 60*33391Sbostic pow_zi(&p1, &a1, b); 6129965Smckusick 62*33391Sbostic p->real = p1.dreal; 63*33391Sbostic p->imag = p1.dimag; 6429965Smckusick } 65*33391Sbostic 66*33391Sbostic #endif /* tahoe */ 67