Commit | Line | Data |
---|---|---|
d6cc74ee DW |
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 | } |