Commit | Line | Data |
---|---|---|
547779a8 WH |
1 | #include "f2c.h" |
2 | ||
3 | #ifdef KR_headers | |
4 | VOID pow_zi(p, a, b) /* p = a**b */ | |
5 | doublecomplex *p, *a; integer *b; | |
6 | #else | |
7 | extern void z_div(doublecomplex*, doublecomplex*, doublecomplex*); | |
8 | void pow_zi(doublecomplex *p, doublecomplex *a, integer *b) /* p = a**b */ | |
9 | #endif | |
10 | { | |
11 | integer n; | |
12 | double t; | |
13 | doublecomplex x; | |
14 | static doublecomplex one = {1.0, 0.0}; | |
15 | ||
16 | n = *b; | |
17 | p->r = 1; | |
18 | p->i = 0; | |
19 | ||
20 | if(n == 0) | |
21 | return; | |
22 | if(n < 0) | |
23 | { | |
24 | n = -n; | |
25 | z_div(&x, &one, a); | |
26 | } | |
27 | else | |
28 | { | |
29 | x.r = a->r; | |
30 | x.i = a->i; | |
31 | } | |
32 | ||
33 | for( ; ; ) | |
34 | { | |
35 | if(n & 01) | |
36 | { | |
37 | t = p->r * x.r - p->i * x.i; | |
38 | p->i = p->r * x.i + p->i * x.r; | |
39 | p->r = t; | |
40 | } | |
41 | if(n >>= 1) | |
42 | { | |
43 | t = x.r * x.r - x.i * x.i; | |
44 | x.i = 2 * x.r * x.i; | |
45 | x.r = t; | |
46 | } | |
47 | else | |
48 | break; | |
49 | } | |
50 | } |