1 | # power function for double precision floating point
|
2 |
|
3 |
|
4 |
|
5 | dpow = ->
|
6 | a = 0.0
|
7 | b = 0.0
|
8 | base = 0.0
|
9 | expo = 0.0
|
10 | result = 0.0
|
11 | theta = 0.0
|
12 |
|
13 | expo = pop_double()
|
14 | base = pop_double()
|
15 |
|
16 | # divide by zero?
|
17 |
|
18 | if (base == 0.0 && expo < 0.0)
|
19 | stop("divide by zero")
|
20 |
|
21 | # nonnegative base or integer power?
|
22 |
|
23 | if (base >= 0.0 || (expo % 1.0) == 0.0)
|
24 | result = Math.pow(base, expo)
|
25 | push_double(result)
|
26 | return
|
27 |
|
28 | result = Math.pow(Math.abs(base), expo)
|
29 |
|
30 | theta = Math.PI * expo
|
31 |
|
32 | # this ensures the real part is 0.0 instead of a tiny fraction
|
33 |
|
34 | if ((expo % 0.5) == 0.0)
|
35 | a = 0.0
|
36 | b = Math.sin(theta)
|
37 | else
|
38 | a = Math.cos(theta)
|
39 | b = Math.sin(theta)
|
40 |
|
41 | push_double(a * result)
|
42 | push_double(b * result)
|
43 | push(imaginaryunit)
|
44 | multiply()
|
45 | add()
|