Commit | Line | Data |
---|---|---|
d6cc74ee | 1 | /* |
a7868ad7 RE |
2 | * Copyright (c) 1980 Regents of the University of California. |
3 | * All rights reserved. The Berkeley software License Agreement | |
4 | * specifies the terms and conditions for redistribution. | |
5 | * | |
6 | * @(#)pow_zi.c 5.1 %G% | |
d6cc74ee DW |
7 | */ |
8 | ||
9 | #include "complex" | |
10 | ||
11 | pow_zi(p, a, b) /* p = a**b */ | |
12 | dcomplex *p, *a; | |
13 | long int *b; | |
14 | { | |
15 | long int n; | |
16 | double t; | |
17 | dcomplex x; | |
18 | ||
19 | n = *b; | |
20 | p->dreal = 1; | |
21 | p->dimag = 0; | |
22 | ||
23 | if(n == 0) | |
24 | return; | |
25 | if(n < 0) | |
26 | { | |
27 | n = -n; | |
28 | z_div(&x, p, a); | |
29 | } | |
30 | else | |
31 | { | |
32 | x.dreal = a->dreal; | |
33 | x.dimag = a->dimag; | |
34 | } | |
35 | ||
36 | for( ; ; ) | |
37 | { | |
38 | if(n & 01) | |
39 | { | |
40 | t = p->dreal * x.dreal - p->dimag * x.dimag; | |
41 | p->dimag = p->dreal * x.dimag + p->dimag * x.dreal; | |
42 | p->dreal = t; | |
43 | } | |
44 | if(n >>= 1) | |
45 | { | |
46 | t = x.dreal * x.dreal - x.dimag * x.dimag; | |
47 | x.dimag = 2 * x.dreal * x.dimag; | |
48 | x.dreal = t; | |
49 | } | |
50 | else | |
51 | break; | |
52 | } | |
53 | } |