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